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CHAPTER  1 


INTRODUCTION 

Finite  dimensional  systems  have  Markovian  or  state  space 
descriptions.  This  is  well  understood.  Similarly,  narrowband  signals, 
such  as  speech,  vibrations,  and  other  natural  phenomena  have  state 
space  descriptions.  In  such  descriptions,  all  relevant  information 
is  carried  in  the  state.  The  Kalman  filtering  problem  is  one  of 
estimating  the  state  of  a  system  based  upon  noisy  measurements  up  to 
time  t.  Kalman  prediction  is  the  process  of  estimating  the  state 
based  upon  noisy  measurements  up  to  time  t-1.  The  Kalman  predictor  is 
actually  a  one  step  predictor. 

In  general  the  Kalman  filter  (predictor)  is  the  best  filter 
among  all  linear  filters.  When  the  noise  processes  are  Gaussian  the 
Kalman  filter  (predictor)  is  the  best  linear  filter  among  all  filters 
linear  or  nonlinear.  The  goodness  is  measured  in  terms  of  the  error 
covariance. 

Among  all  state  space  models  of  a  process  there  is  one  type  of 
model  which  is  of  special  interest.  This  is  the  so-called  Innovations 
representation.  The  name  comes  from  the  fact  that  the  model  is 
driven  by  the  innovations  process  of  the  output.  It  so  happens  that 
the  innovations  process  of  the  Kalman  predictor  (filter)  is  the  same 
as  the  innovations  process  of  the  signal  model  when  there  is  no  mea¬ 
surement  noise.  The  Kalman  predictor  (filter)  is  an  exact  inverse  of 
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this  kind  of  signal  model.  This  fact  will  be  used  later  to  define  the 
Kalman  predictor  (filter). 

In  this  thesis  we  restrict  our  attention  to  autoregressive 
moving  average  (ARMA)  processes.  We  study  the  Kalman  filter/ 
predictor  both  from  a  theoretical  and  from  a  practical  point  of  view. 
Numerical  results  are  presented  for  the  Kalman  filter  implemented  in 
floating  point  and  fixed  point  arithmetic.  What  follows  is  a  short 
overview  of  each  chapter. 

Chapter  2  is  devoted  to  signal  models  for  autoregressive  moving 
average  processes  and  their  properties.  The  first  two  signal  models 
we  derive  are  the  Markovian  representations  #1  and  #2.  They  are 
linear  time-invariant  state  space  models.  Innovations  representa¬ 
tions  #1  and  #2  are  on  the  other  hand  linear  time-varying  state  space 
models.  The  state  vectors  of  these  models  are  vectors  of  s-step 
predictors  based  on  the  infinite  past  for  the  Markovian  state  space 
models  and  based  on  the  finite  past  for  the  Innovations  representa¬ 
tions. 

In  Chapter  3  we  introduce  the  Kalman  predictor  (filter)  as  a 
sort  of  inverse  to  its  assumed  signal  model.  It  is  the  Innovations 
representation  H  for  the  Kalman  predictor  and  the  Innovations  repre¬ 
sentation  #2  for  the  Kalman  filter.  Also  we  look  at  some  of  the 
properties  of  the  Kalman  predictor  and  filter.  The  Kalman  gain  vec¬ 
tors  are  composed  of  the  first  n+1  values  of  a  time-varying  unit 
pulse  response.  A  so-called  Impulse  Response  Algorithm  is  derived 
which  allows  us  to  calculate  the  unit  pulse  response  recursively  and 
obtain  the  Kalman  gain  vectors.  Using  the  fact  that  the  underlying 
process  is  an  ARMA  (p,q)  a  substantial  savings  in  calculation  is 
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achieved.  The  result  Is  the  so-called  Fast  Kalman  Algorithm. 

Besides  being  fast,  this  algorithm  is  very  well  suited  for  fixed 
point  implementation.  It  is  shown  that  all  the  elements  of  the  gain 
vectors  are  less  than  or  equal  to  one  in  magnitude. 

Chapter  4  is  devoted  to  some  of  the  problems  which  must  be 
addressed  when  implementing  the  Kalman  predictor  and  filter  using 
fixed  point  arithmetic.  Scaling  of  variables  is  used  to  limit  the 
probability  of  overflow.  This  is  seldom  a  problem  in  floating  point 
arithmetic,  but  in  fixed  point  arithmetic  it  is  nearly  always  a  prob¬ 
lem.  The  effects  of  roundoff  can  be  modeled  as  additive  white  noise 
sources  and  incorporated  into  the  governing  equations  for  the  Kalman 
predictor  and  Kalman  filter.  Other  effects  as  a  result  of  coeffi¬ 
cient  truncation  and  input  quantization  are  also  important  but  not 
considered  in  this  thesis.  We  do  analyze  computational  complexity  in 
terms  of  number  of  arithmetic  operations. 

In  Chapter  5  numerical  results  are  presented.  The  Kalman 
filter  is  implemented  using  both  floating  point  and  fixed  point 
arithmetic  and  the  results  are  compared.  As  it  turns  out  the  fixed 
point  results  compare  very  favorably  with  the  floating  point  results 
even  though  the  most  simple  scaling  procedure  is  used.  As  mentioned 
before,  the  Fast  Kalman  Algorithm  is  very  we’l  suited  for  fixed  point 
implementation  and  the  numerical  results  confirm  this.  In  Chapter  6 
we  summarize  the  main  results  of  this  thesis  and  suggest  further 
research. 


CHAPTER  2 


SIGNAL  MODELS  FOR  ARMA  PROCESSES 

This  chapter  deals  with  state  space  representations  of  an 
autoregressive  moving  average  (ARMA)  process.  Markovian  representa¬ 
tions  are  linear  time-invariant  state  space  models.  Their  correspond¬ 
ing  state  vectors  are  vectors  of  s-step  ahead  predictors  based  on  the 
infinite  past.  Innovations  representations  are  a  linear  time  varying 
state  space  models  with  state  vectors  of  s-step  ahead  predictors 
based  on  the  finite  past. 

In  section  2.1  we  derive  two  closely  related  Markovian  state 
space  models  which  we  call  Markovian  representations  #1  and  #2  (MRl 
and  MR2).  Section  2.2  Is  devoted  to  some  of  the  properties  of  the 
Markovian  representations.  Then  in  section  2.3  two  Innovations  rep¬ 
resentations  (IRl  and  IR2)  are  Introduced  as  time-varying  Markovian 
counterparts.  In  section  2.4  the  properties  of  the  Innovations  rep¬ 
resentations  are  studied  and  their  properties  are  compared  with  those 
of  the  Markovian  state  space  models. 

2.1  Markovian  Representations 

In  this  section  we  derive  two  closely  related  state  space 
representations  of  an  autoregressive  moving  average  process. 

Autoregressive  moving  average  process.  A  second  order  sta¬ 
tionary  (SOS)  time  series  y(t)  is  called  an  autoregressive  moving 
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average  process  of  order  (p»q)»  (ARMA(p,q)),  if  it  satisfies  the 
equation 

y(t)  +  a(l)y(t  -  1)  +  .  .  .  +  a(p)y(t  -  p)  = 

b(0)u(t)  +  b(l)u(t  -  1)  +  .  .  .  +  b{q)u(t  -  q)  (2.1.1) 
In  operator  form  this  equation  is 

A(z)  y(t)  *  B(z)  u(t)  (2.1.2) 

where 

A(z)  =  1  +  a(l)z‘^  +  .  .  .  +  a(p)z’P  (2.1.3) 

B(z)  =  b(0)  +  b(l)z"^  +  .  .  .  +  b(q)z'‘’  (2.1.4) 

The  input  u(t)  is  a  white  noise  process: 

(a  2  i=0 

ECu(t)]  =  0  and  E[u(t)u(t  +  i)]  =  \  ” 

(  0  i  not  zero 

The  zeros  of  A(z)  and  B(z)  are  assumed  to  lie  inside  the  unit  circle. 
See  Figure  2.1.1  for  a  block  diagram.  The  unit  pulse  response  of  the 
ARMA(p,q)  process  is 

h(t)  +  a(l)h(t  -  1)  +  .  .  .  +  a(p)h(t  -  p)  = 

b(0)fi(t)  +  b(l)6(t  -  1)  +  .  .  .  +b(q)6(t  -  q)  (2.1.5) 
Best  linear  predictor.  The  fact  the  zeros  of  A(z)  and  B(z)  lie 
inside  the  unit  circle  allows  us  to  write 

00 

y(t)  =  A"^(z)B(z)u(t)  =  E  h(i)u(t  -  1)  = 

i=0 

y(l/t~l)  +  h(0)u(t)  (2.1.6) 

and 

CO 

u(t)  =  B‘^(z)A(z)y(t)  =  E  d(i)y(t  -  i )= 

i=0 

d(0)[y(t)  -  y(l/t-l)l  (2.1.7) 

where  y(l/t-l)  is  the  best  (one  step)  linear  '>redictor  i  UP)  for 
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y(t)  given  y(t-i)  i  =  1,2,  ....  It  is  defined  as 

d(0)y(l/t-l)  =  -d(l)y(t-l)  -  d(2)y(t-2)  -  .  .  .  = 

00 

-  i:d(i)y(t-i) 
i=l 

where  the  coefficients  d(i)  i=l,2,3,  .  .  .  minimize 
E[y(t)  -  y(l/t-l)]2 

In  terms  of  u(t-i)  i=l,2,3,  .  .  .  the  BLP  is 
00 

y(l/t-l)  =  Z  h(i)u(t-i) 
i=l 

and  the  prediction  error  e(t)  is  then 

e(t)  =  y(t)  -  y(l/t-l)  =  h(0)u(t) 

More  generally  the  best  linear  s-step  predictor  given  y(t-i ) 
i =1,2,3,  .  .  .  is 


y{s/t-l)  =  Z  h(s+i )u{t-l-i )  s=  1,2,3,  .  .  . 
i=0 

or 


y(s/t)  =  Z  h(s+i)u(t-i)  s=  1,2,3,  .  .  . 
i=0 

Markovian  Representation  ^1.  We  can  now  write 
y(s/t)  =  h(s)u(t)  +  h(s+l)u(t-l)  +  .  .  . 
y(s+l/t-l)  =  h(s+l)u(t-l)  +  .  .  . 

and  after  subtracting  the  second  equation  from  the  first 
y(s/t)  =  y(s+l/t-l)  +  h(s)u(t) 

Using  the  fact  that  (see  Box  and  Jenkins  [31]) 

y(s/t)  +  a(l)y(s-l/t)  +  .  .  .  +  a(n)y(s-n/t)  =  0  s>=n 


(2.1.8) 

(2.1.9) 

(2.1.10) 

(2.1.11) 

(2.1.12) 

(2.1.13) 

(2.1.14) 


y(n/t-l)  =  -a(l)y(n-l/t-l)-  .  .  .  -a(n)y(l/t-l ) 


(2.1.15) 
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where  n  =  max(p,q)  we  can  write  down  a  state  space  representation 
using  equation  (2.1.13)  for  s=l,2,  .  .  .,  n  and  (2.1.15) 


x(t+l/t)  =  Ax(t/t-l)  +  h(l)u(t)  (2.1.16) 

y(t)  =  y(l/t-l)  +  h(0)u(t) 

=  c'x(t/t-l)  +  h(0)u(t)  (2.1.17) 

Here 

x'(t/t-l)  =  [y(l/t-l).y(2/t-l),  .  .  .  y(n/t-l)]  (2.1.18) 

h’(l)  =  [h(l),h(2),  .  .  .,h(n)]  (2.1.19) 

c‘  =  (1  0  0  ...  0)  (2.1.20) 


The  output  y(t)  is  simply  the  predicted  value  plus  the  prediction 
error.  See  Figure  2.1.2  for  a  block  diagram.  We  will  refer  to  equa¬ 
tions  (2.1.16)  and  (2.1.17)  as  the  Markovian  representation  #1.  Note: 
The  state  vector  consists  of  s-step  ahead  predictors  for 


s=l,2,3,  .  .  .  n  based  on  the  infinite  past: 
x'(t/t-l)  =  Cy(l/t-l),y(2/t-l),  ....  y(n/t-l)] 


Markovian  Representation  #2.  Now  we  write 


y(s/t+l)  =  h(s)u(t+l)  +  h(s+l)u(t)  +  .  .  . 


y(s+l/t)  =  h(s+l)u(t)  +  .  .  . 


and  after  subtracting  the  second  equation  from  the  first; 
y(s/t+l)  =  y(s+l/t)  +  h(s)u(t+l) 


(2.1.22) 


Fit.  2.1.2.  Markovian  representation  #1. 


Note  that 

y(n/t)  =  -a(n)y(n-l/t)-  .  .  .  -a(l)y(0/t) 

We  can  write  for  s=0,l,  .  .  .,n-l  another  state  space 

representation 

)((t+l/t+l)  =  A2<(t/t)  +  ]i(0)u(t+l) 

(2.1.23) 

y(t)  =  c'x(t/t) 

(2.1.24) 

Here 

x'(t/t)  =  [y(0/t),y(l/t),  ....  y(n-l/t)] 

(2.1.25) 

ll'(O)  =  [h(0).h(l),  .  .  .,  h(n-l)] 

(2.1.26) 

and  c  and  A  are  defined  as  before.  The  output  y(t)  is  the  first 
element  in  the  state  vector  since  y(0/t)  =  y(t)  by  definition.  We 
will  refer  to  equations  (2.1.23)  and  (2.1.24)  as  the  Markovian  repre¬ 
sentation  #2.  This  representation  is  the  same  as  derived  in  the  paper 
by  Akaike  [2]  and  the  previous  derivation  follows  his  closely. 

Figure  2.1.3  shows  a  block  diagram  of  the  Markovian  representation  #2. 
Note: 

The  state  vector  consists  of  s-step  ahead  predictors  for 
s=0,l,2,  .  .  .,  n-1  based  on  the  infinite  past: 
x'(t/t)  =  [y(0/t),y(l/t),  ....  y(n-l/t)] 

The  relation  between  the  two  state  vectors  x(t/t-l)  and  x(t/t)  is 
simply 

x(t/t-l)  =  A  x(t-l/t-l)  (2.1.27) 

Also  note  that 


h(l)  =  A  h(0)  +  d  b(n) 
d*  =  [0  0  ...  1] 
which  follows  from  equation  (2.1.5). 


(2.1.28) 

(2.1.29) 
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2.2  Properties  of  the  Markovian  Representations 


First  and  second  order  statistics  of  the  two  Markovian 
representations  are  the  main  subject  of  this  section.  Let's  begin 
with  the  Markovian  representation  #1. 


Markovian  Representation  H 

Expected  value  of  the  state  vector.  The  MRl  is  described  by 


x(t+l/t)  =  Ax(t/t-l)  +Jh(l)u(t) 

y(t)  =  c'x(t/t-l)  +  h(0)u(t) 

Rewrite  equation  (2.2.1)  as 

t  ^  i  1 

x(t/t-l)  =  a’^x(0/-1)  +  E  A^  ^h(l)u(t-i) 

i=l 


(2.2.1 

(2.2.2 


=  Z  A^"^  h(l)u(t-i) 


(2.2.3 


Taking  the  expectation  on  both  sides  of  equation  (2.2.3) 


E[x(t/t-l)]  =  Z  a’"^  h(l)E[u(t-i)]  =  0 
i=l 


(2.2.4 


since 


E[u(t-i )]  =  0  V  i  >  =  0 


Unit  pulse  response, 


L  The  first  n+1  values  of  the  unit 


pulse  response  are  determined  from  equation  (2.1.5): 


h(t)  =  b(t)  -  Z  a(i)h(t-i)  t=0,l,  .  .  .,  n 
i=l 

where  n=max(p,q)  as  before.  For  t  >  n,  h(t)  is  simply 


(2.2.5 


h(t)  =  -  Z  a(i )h(t-i )  t  >  n 


(2.2.6 


t>=l 


I  c'A’''^h(l) 

h'(l)  =  [h(l).h(2).  .  .  ..  h(n)] 

This  result  can  also  be  obtained  by  writing 
y(t)  =  c'j<(t/t-l)  +  h(0)u(t) 
and  substitute  equation  (3.2.3)  for  x(t/t-l): 

y(t)  =  h(0)u(t)  +  .^j£‘A’‘^Jl(l)u(t-i) 

State  variance.  Q(t).  The  variance  of  state  vector  is  given  by 
Q(t+1)  =  E[x(t+l/t)x'(t+l/t)] 

E[(Ax(t/t-l)+h(l)u(t))(Ax(t/t-l)+h(l)u(t))'] 

-  AQ(t)A'  +  h(l)ll'(l) 

Given  initial  condition  in  the  infinite  past,  Q(-«»),  it  must  hold  that 
Q(0)  =  A*Q(0)A'  +  a^^2h(l)h‘(l)  (2.2.9) 

Q(t)  =  (5(0)  for  t>=0 

State  covariance,  R(t).  The  state  covariance  is  given  by 

R(t)  =  [x(s/s-l)x'(s+t/sn-l)] 

t  t  i  1 

E[x(s/s-l)(A^x(s/s-l)  +  I  A^  ^  h(l)u(t+s-i))'] 

i=l 

=  E[x(s/s-l)x.' (s/s-l)](A^)' 

=  A^Q(s) 

=  A^Q(O)  (2.2.10) 

Output  covariance,  r(t).  The  covariance  of  y(t)  is  determined 
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r(t)  =  E[y(s)y'(s+t)] 

E[(£'x(s/s-1)  +  h(0)u(s))(£'x(s+t/s+t-l)  +  h(0)u(s+t))'] 
=  £'E[)c(s/s-1)_x'  (s+t/s+t-l)]£+h(0)E[u(s)x'  (s+t/s+t-l)]£ 


£’R(t)  £  +  0^2  h(0)h{t) 


(2.2.11) 


Here  we  have  used 


E[u(s)x'(s+t/s+t-l)]  =0^2  (A^‘V(1))' 

=  r.' \ 


(2.2.12) 


h(t)  =  £'A‘'‘^h(l) 

Also  note  that  r(t)  =  r(-t). 

By  using  similar  techniques  we  derive  the  corresponding 
properties  for  the  Markovian  representation  #2. 


Markovian  Representation  #2 

Expected  value  of  the  state  vector.  The  MR2  is  described  by 
)((t+l/t+l)  =  Ax(t/t)  +  h^(0)u(t+l) 
y(t)  =  £'x(t/t) 


Writing 


t  . 


x(t/t)  =  A^'^^x(-1/-1)  +  E  A’h(O)  u(t-i) 

i=0 


=  E  A^fi(O)  u(t-i ) 
i=0 

and  taking  the  expectation  on  both  sides: 


(2.2.13) 


E[2i(t/t)]  =  E[  E  A’h^(0)u(t-i )]  =  0 
i=0 


(2.2.14) 


The  unit  pulse  response.  h(t).  This  follows  by  inspection 
from  equation  (2.2.13)  and  is  given  by 


h(t)  = 


(2.2.15) 


c'A^  h(0)  t>0 
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State  variance.  g(t).  The  state  variance  is  given  by 


$(t+l)  =  AC(t)A'  +  ay2h(0)h'(0)  (2.2.16) 

The  stationary  variance  satisfies 

e(0)  =  Ae(0)A'  +  o^^2h(0)h'(0)  (2.2.17) 

State  covariance.  i?(t).  The  state  covariance  is  given  by 
i?(t)  =  A^$(0)  t>=0  (2.2.18) 

Output  covariance.  r(t).  The  output  covariance  of  MR2  equals 
that  of  the  MRl  and  is  given  by 

r(t)  =  c'ff(t)  c  t>=0  (2.2.19) 


Table  1  summarizes  the  basic  properties  of  the  two  Markovian  represen¬ 
tations  we  have  derived  is  this  section. 

2.3  Innovations  Representations 

In  the  time  invariant  Markovian  state  space  representations 
studied  in  the  previous  sections  the  non-zero  Initial  conditions  are 
random  vectors: 

x(0/-l)  :  N(0,Q(0)) 
x(0/0)  :  N(0,e(0)) 

where  Q(0)  and  Q{0)  are  determined  from 
Q(0)  =  AQ(0)A'  +  ay2h(l)h'(l) 

5(0)  =  Ag(0)A'  +  o^^2h(0)]l'(0) 

One  way  of  deriving  the  innovations  representations  for  an  ARMA(p,q) 
process  is  to  introduce  them  as  inverses  of  a  Kalman  predictor/filter. 
This  is  done  by  Anderson  and  Moore  [4].  It  may  be  done  more  directly 
by  treating  them  as  time  varying  counterparts  of  the  Markovian  repre¬ 
sentations  as  done  by  Gueguen  and  Scharf  [1].  The  Kalman  predictor/ 
filter  is  introduced  as  a  kind  of  inverse  of  the  Innovations  represen¬ 
tations.  I  choose  the  latter  way. 


■  h{n)] 


Replacing  ji(l)  by  k(t),  ii(0)  by  fe(t+l)  and  by  v(t)  in  the 
Markovian  representations  we  get  time-varying  representations  called 
the  Innovarions  representations. 

Innovations  representation  #1.  The  defining  equations  are 
x(t+l/t)  =  Ax(t/t-l)  +  k(t)u(t)  (2.3.1) 


x(t+l/t)  =  Ax(t/t-l)  +  jc(t)u(t) 
y(t)  =  c'x(t/t-l)  +  h(0)u(t) 
where  Jc(t)  is  a  time-varying  vector 

k'(t)  =  Ck^l).  k^(2).  ....  k^(n)] 
and  u(t)  is  the  innovations  sequence 
u(t)  :  N(0,v(t)) 


(2.3.2) 


(2.3.3) 


(2.3.4) 


Innovations  representation  #2.  The  defining  equations  are 


x(t+l/t+l)  =  Ax(t/t)  +  fe(t+l)u(t+l) 
y(t)  =  c'x(t/t) 

where  kit)  is  a  time  varying  vector 

k'M  •  U‘(l),  k*(2) . *:*(")] 

and  u(t)  the  innovations  sequence 
u(t)  :  N(0,v(t)) 


(2.3.5) 

(2.3.6) 


(2.3.7) 


(2.3.8) 


The  initial  conditions  on  the  state  vectors  )^(t/t-l)  and  x(t/t)  are 
chosen  to  be  identically  zero: 

x(0/-l)  =  0  (2.3.9) 

x(-l/-l)=0  (2.3.10) 

The  definitions  for  A  and  c  are  the  same  as  for  the  Markovian 
representations. 

The  Innovations  representations  are  the  time-varying  state 
space  models  that  start  from  zero  initial  conditions  rather  than  from 
random  initial  conditions  as  in  the  Markovian  representations  [1]. 

The  state  vectors  of  the  Innovations  representations  are  still 
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vectors  of  s-step  predictors  but  based  on  the  finite  past: 

x'(t/t-lKy(l/t-l).y(2/t-l).  .  .  ..  y(n/t-l)]  (2.3.11) 

x.'(t/t)  -  [y(0/t).y(l/t) . y{n-l/t)]  (2.3.12) 

Figures  2.3.1  and  2.3.2  show  block  diagrams  of  the  two  Innovations 
representations. 

In  section  2.4  we  show  how  to  choose  Jc(t),  k(t)  and  v(t)  in 
order  to  generate  stationary  data  y(t)  with  the  same  correlations 
sequence  as  that  of  the  Markovian  representations. 

2.4  Properties  of  the  Innovations  Representations 

As  in  section  2.2  we  derive  some  first  and  second  order  statis¬ 
tics  for  the  two  Innovations  representations  and  then  show  what 
conditions  must  be  satisfied  to  match  the  output  covariance  of  the 
Innovations  representations  to  the  one  of  the  Markovian  representa¬ 
tions.  We  start  with  the  Innovations  representation  #1. 

Innovations  Represe  M  wa  tion  #1 

Expected  value  of  the  state  vector.  The  IRl  is  described  by 
£(t+l/t)  =  A)((t/t-l)  +  jc(t)u(t) 
y(t)  *  c'x(t/t-l)  +  h(0)u(t) 

The  first  equation  can  be  written  as 

t  ^  i-1 

x(t/t-l)  =  A^x(0/-1)  +  r  a’  ^k(t-i)u(t-i) 

i»l 


CO 

«  E  A^"^k(t-i)u(t-i)  (2.4.1) 

i=l 

Taking  the  expectation  on  both  sides  we  get 

E[x(t/t-l)]  =  Z  A^"^k(t-i)E[u(t-i)]  =  0 
i=l 


(2.4.2) 


19 


a 
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Unit  pulse  response.  We  know  that  the  unit  pulse  re¬ 

sponse  must  satisfy 


y(t)  =  Z  h^(i)u(t-i) 
i=0 


(2.4.3) 


Now 


y(t)  =  c'x{t/t-l)  +  h(0)u{t)  (2.4.4) 

Substituting  equation  (2.4.1)  for  )((t/t-l)  in  (2.4.4)  to  obtain 

00 

y(t)  =  h(0)u(t)  +  Z  A^"^k(t-i)u(t-i)  (2.4.5) 

i=l 

Compare  equations  (2.4.3)  and  (2.4.5)  to  obtain  the  unit  pulse 
response: 

0  i<0 


h*(1)  ' 


h(0) 


i=0 


(2.4.6) 


c'A^’^k(t)  i>0 


Note  that  h^(i)  is  the  response  at  time  t+i  to  an  impulse  applied  at 
time  t. 

State  variance.  Q(t).  The  variance  of  the  state  vector  is 
given  by 

Q(t+1)  *  E[x(t+l/t)x'(t+l/t)]  = 

E[(A2<(t/t-l)+k(t)u(t))(A2((t/t-l)+Jc(t)u(t)) '] 

=  AQ(t)A'  +  v(t)k(t)k'(t)  (2.4.7) 

Here 

Q(0)  =  0  (2.4.8) 

v(t)  =  E[u(t)u'(t)]  (2.4.9) 


state  covariance.  R^(k).  The  state  covariance  is  given  oy 
R^(k)  =  [x{t/t-l)x'(t+k/t+k-l)] 

=  E[x(t/t-l)(A  x(t/t-l)  +  E  ^k(t+k-i)u(t+k-i))'] 

i=l 

=  A'^Q(t)  (2.4.10) 

since 

E[)<(t/t-l)u' (t+k-i )]  *  0  for  i=l,2,  .  .  .,  k 

Output  covariance,  r^(k).  The  covariance  of  y(t)  is  detentiined 

from 

r^(k)  =  E[y(t)y'(t+k)] 

=  E[(c')<(t/t-l  )+h(0)u{t) )  (£'><(t+k/t+k-l  )+h(0)u(t+k) ) '  ] 

=  c'R^(k)  c  +  v(t)h(0)h^(k)  (2.4.11) 

This  follows  from 

E[u(t)x'(t+k/t+k-l)]  =  v(t)(A'^’^k(t))'  (2.4.12) 

and  equation  (2.4.6). 

Innovations  Representation  #2 

Expected  value  of  the  state  vector.  The  IR2  is  described  by 
2i(t+l/t+l)  =  Ax(t/t)  +  /c(t+l)u(t+l) 
y(t)  =  c'x(t/t) 

The  state  equation  can  be  written  as 

x(t/t)  =  A^'*’^x(-1/-1)  +  I  A’fe(t-i)u(t-i)  (2.4.13) 

i=0 

Taking  the  expectation  on  both  sides  results  in 

t  . 

E[x(t/t)]=E[A^'^^x(-l/-l)  +  ^  A’^(t-i)u(t-i)]=0 

i=0 


(2.4.14) 
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Unit  pulse  response,  h  (i).  We  have 

y(t)  =  c'2((t/t)  =  c'  Z  A^fe(t-i  )u(t“i ) 

i=0 

so  the  unit  pulse  response  is  simply 


h^i)  = 


(2.4.15) 


c'A  k(t) 


By  comparison  with  equation  (2.4.6)  the  relation  between  J<(t)  and 
fe(t)  is 


k(t)  =  A  fe(t) 


(2.4.16) 


State  variance,  o(t).  The  state  variance  is  determined  from 
«(t)=E[x(t/t)x'(t/t)] 

=  E(Axit-l/t-l)+fe(t)u(t))(Ax(t-l/t-l)+fe(t)u(t))'] 

=  A$(t-1)A'  +  v(t)fe(t)fe‘(t)  (2.4.17a) 


^(t)  =  Q(t)  +  v(t)/c(t)k'(t) 

Q(t)  is  the  state  variance  of  the  Innovations  rep.  #1. 


(2.4.17b) 


State  covariance,  E  (k).  The  state  covariance  is  given  by 
/?^(k)=E[x(t/t)x' (t+k/t+k)]  = 


E[x(t/t)(A'^x(t/t)  +  Z  A^fc(t+k-i)u(t+k-i))'] 

i=0 

=  A^(t) 


(2.4.18) 


since 


E[x(t/t)u’(t+k-i )]  =  0  for  i=0,l,2 . k-1 

Output  covariance,  r^(k).  The  covariance  of  y(t)  is  detemiined 


I 
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r^(k)  =  E[y(t)y'(t+k)]  =  c'E[x(t/t)x' (t+k/t+k)]c 

=  c'  R^{k)  c  (2.4.19) 

Table  2  summarizes  the  basic  properties  of  the  two  Innovations  repre¬ 
sentations  we  have  derived  in  this  section. 

The  output  covariance  r^(k)  of  an  Innovations  representation  is 
generally  non-stationary  but  we  would  like  to  make  it  stationary  or 
more  specifically  match  it  to  that  of  the  corresponding  Markovian 
representation.  For  the  Innovation  representations  #1 
r^(k)  =  c'R^(k)c  +  v(t)h(0)h^{k) 

=  c'A'^Q(t)c  +  h(0)c‘A'^'^k(t)v(t)  (2.4.20) 

Gueguen  and  Scharf  [1]  show  that  by  choosing 

k(t)v(t)  =  -AQ(t)c  +  AQ(0)c  +  o^j2h(l)  (2.4.21) 

Q(0)  =  AQ(0)A'  +  a^^2h(i)h'(i) 

r^(k)  is  made  equal  to  r(k)  for  the  Markovian  repre' entation  #1. 

Plug  equation  (2.4.21)  into  (2.4.20) 

r'^(k)=h(0)c'A'^Q(0)c+a^j2h(0)c'A'^'^h(l)+(l-h(0))c'A'^Q(t)c 

or 

r^(k)  =  c'a'^Q(0)c  +  a^^2h(k)  =  r(k)  (2.4.22) 

The  variance  of  the  innovations  sequence  u(t),  v(t),  becomes 

v(t)=E[u(t)u'(t)]=E[(y(t)-c'x(t/t-l))(y(t)-c'x(t/t-l))'] 

=  E[y(t)y' (t)]-c'E[x{t/t-l)x' (t/t-l)]c 
=  r(0)  -  c'Q(t)c 

=  c'(0(O)-Q(t))c  +  (2.4.23) 

since  r(0)  =  c'Q(0)c  + 

The  output  covariance  for  Innovations  representations  #2  is 
r^(k)=c'A'^$(t)c=c'A^A5(t-l)A'+v(t)k(t)^'(t))c 
=  c'A'^'^^(3(t-l)A'c  +  c'h(0)A'^fe(t)v(t) 


(2.4.24) 


Table  2.  Governing  equations  for  the  Innovations  representation  #1  and  HZ. 


(l)e--  •  •  •  (u)e- 


(2.4.25) 


Now  we  choose  fe(t)v(t)  as  follows 

fe(t)v(t)  =  -AQ(t-l)A'c  +  S(0)c 
where  5(0)  satisfies 

5(0)  =  A*5(0)*A'  +  Oy2h(0)h'(0) 

Plug  equation  (2.4.25)  into  (2.4.24) 

r^(k)  =  c'A‘^‘^^5(t-l)A’c+c'h(0)A’^[-A5(t-l)A'+5(0)]c 

=  c’h(0)A'^5(0)c  +  (l-h(0))c'A'^‘^^5(t-l)A’c  (2.4.26) 

or 

r^(k)  =  r(k)  =  c'a'^5(0)c  (2.4.27) 

The  variance  of  the  innovations  sequence  u(t),  v(t),  is  now 
v(t)  =  E[u(t)u'(t)]  =  r(0)  -  c'A5(t)A'c 

=  c'A[5(0)-5(t-l)]A'c  +  (2.4.28) 

The  main  point  is  this.  It  is  possible  to  generate  stationary  data, 
using  the  Innovations  representations,  having  the  same  first  and 
second  order  statistics  as  the  Markovian  representations.  For  the 
Innovation  representation  #1  this  is  done  by  choosing 
k(t)v(t)  =  -A*Q(t)c  +  A*Q(0)c  +  o^j2h(l) 
h(0)  =  1 

and  for  the  Innovations  representation  #2  by  choosing 
k(t)v(t)  =  -A*5(t-l)*A'c  +  5(0)c 
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CHAPTER  3 


FAST  KALMAN  FILTERING  FOR  ARMA  PROCESSES 

I 

Using  the  two  time-varying  signal  models  introduced  in  Chapter 
2  we  derive  the  Kalman  predictor  and  filter  as  inverses  of  these  sig- 
I  nal  models.  In  the  case  of  the  Kalman  predictor  the  assumed  signal 

model  is  the  Innovations  representation  #1  and  for  the  Kalman  filter  it 
is  the  Innovations  representation  #2.  Anderson  and  Moore  [2]  intro¬ 
duced  the  Innovation  representation  #1  as  a  sort  of  inverse  to  a  given 
Kalman  predictor  but  the  idea  of  doing  the  opposite  is  attributed  to 
Gueguen  and  Scharf  [1].  They  derived  the  Innovations  representation 
#1  as  the  time-varying  counterpart  of  the  Markovian  representation  #1. 
In  section  3.1  and  3.2  we  introduce  the  Kalman  predictor  and  Kalman 
filter,  respectively. 

The  calculations  of  the  Kalman  gain  vectors,  k(t)  and  fc(t),  are 
the  main  computational  problems  in  implementing  the  Kalman  predictor 
and  filter.  Solving  a  Ricatti  type  of  equation  is  not  the  most  effi¬ 
cient  way  and  a  better  algorithm  is  desired.  The  problem  of  solving 
the  Kalman  gain  vectors  translates  into  that  of  calculating  the  first 
n+1  values  of  the  time-varying  impulse  response  for  eact  t.  A  so- 
called  Impulse  Response  Algorithm  can  be  used  to  calculate  the  impulse 
response  recursively.  The  number  of  steps  must  be  known  ahead  of  time. 
This  drawback  is  bypassed  by  using  the  fact  that  the  underlying  pro¬ 
cess  is  an  ARMA(p,q).  This  gives  the  so-called  Fast  Kalman  Algorithm 
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which  calculates  the  Kalman  gain  vectors  recursively.  Dugre'  [7] 
derived  this  algorithm  for  the  predictor  gain  vector  k(t)  and  we 
extend  it  for  the  filter  gain  vector  also.  This  is  the  main  objective 
of  section  3.3. 

Properties  of  the  Kalman  predictor  and  filter  are  the  subject 
of  section  3.4.  Both  the  signal  models  and  their  corresponding  pre¬ 
dictor  and  filter  structures  are  asymptotically  time  invariant  and  in 
the  case  of  the  Innovations  representations  #1  and  #2  they  converge  to 
the  time  invariant  Markovian  representations  #1  and  #2,  respectively. 
The  Kalman  predictor  and  filter  converge  to  their  corresponding  Mark¬ 
ovian  predictor  and  filter  structures. 

3.1  The  Kalman  Predictor 

The  (one  step)  prediction  problem  is  one  of  estimating  the 
state  vector,  x(t),  of  the  signal  model  given  measurement  data  up  till 
time  t-1.  For  more  information  on  the  general  problem  of  optimal  pre¬ 
diction  see  references  [2],  [3],  [4],  [5],  and  [6].  We  consider  two 
different  cases: 

1)  Noise-Free  Measurement  Process:  z{t)  =  y(t) 

2)  Noisy  Measurement  Process:  z(t)  =  y(t)  +  n(t) 

Here  n(t)  is  a  sequence  of  independent  N(0,o^^)  random  variables. 

Noise-free  measurement  processes.  The  Innovations  rep.  #1  is 
described  by 

x(t+I/t)  =  Ax(t/t-l)  +  ]<{t)u(t);  x(0/-l)=0  (3.1.1) 

y(t)  -  c'x(t/t-l)  ^  u(t);  h(0)-l  (3.1.2) 

Solve  for  u(t)  in  equation  (3.1.2): 

u(t)=y(t)-c'x(t/t-l)=y(t)-y(t)=z(t)-y(t) 


(3.1.3) 
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Plug  u(t)  into  equation  (3.1.1)  to  obtain  the  Kalman  predictor: 

x(t+l/t)=Ax(t/t-l)+k(t)[z{t)-y(t)];  x(0/-l)=0  (3.1.4) 

y(t)  =  c'x(t/t-l)  (3.1.5) 

A  derivation  similar  to  this  one  can  be  found  in  [1].  In  [2]  the 
Innovations  representation  #1  is  derived  as  an  inverse  to  the  Kalman 
predictor.  Figure  3.1.1  shows  the  Innovations  representation  #1  (the 
signal  model)  with  the  Kalman  predictor.  In  the  case  of  a  noise-free 
measurement  process  the  Kalman  predictor  reconstructs  the  innovations 
sequence  of  the  signal  model  (the  Innovation  rep.  #1)  and  estimates  the 
state  vector  exactly.  Remember,  though,  that  the  perfectly  estimated 
state  consists  of  imperfect  estimate  of  the  underlying  sequence  y(t). 
The  Kalman  predictor  can  also  be  said  to  be  a  whitening  filter  if  we 
think  of  u(t)  as  the  output.  The  input  is  a  correlated  data  sequence 
z(t)  and  the  output  is  uncorrelated  data  sequence  u(t).  The  defining 
equations  for  the  Kalman  predictor  are  then: 


x(t+l/t)=Ax(t/t-l)+j^(t)[z(t)-y(t)];  x(0/-l)=0 

y(t)  =  c'x(t/t-l) 

k(t)v(t)  =  A[Q(0)-Q(t)]c'  +  a^^2h(l) 

(3.1.6) 

v(t)  =  r(0)  -  c'Q(t)c 

(3.1.7) 

Q(t+1)  =  AQ(t)A'  +  v(t)k(t)k'(t);  Q(0)=0 

(3.1.8) 

Q(0)  =  AQ(0)A'  +  a^2h(l)h'(l) 

(3.1.9) 

r(0)  =  c'Q(0)c  + 

(3.1.10) 

These  are  basically  the  equations  derived  in  section  2.4. 

Noisy  measurement  process.  When  there  is  additive  noise,  the 


input  to  the  Kalman  predictor  is 
2(t)  =  y(t)  +  n(t) 


(3.1.11) 
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where  n(t)  is  a  sequence  of  independent  N(0,a^*)  random  variables. 
The  noise  added  to  the  output  is  only  going  to  affect  the  output  co- 
variance  of  the  signal  model.  The  output  covariance  becomes 
r(k)  =  E[z(t)2'(t+k)]  *  E[(y(t)+n(t))(y(t)+n(t))'] 


=  r(k)  +  a^H{k)  (3.1.12) 

and  for  k=0, 

r(0)  =  r(0)  +  0^2  (3.1.13) 

Then  the  defining  equations  for  the  Kalman  predictor  are 

x(t+l/t)  =  Ax(t/t-l)  +  £(t)u(t);  x(0/-l)=0  (3.1.14) 

G(t)  =  2(t)  -  y(t)  (3.1.15) 

y(t)  =  c'£(t/t-l)  (3.1.17) 

k(t)v(t)  =  A[Q(0)-Q(t)]c'  +  a^^^hd)  (3.1.18) 

v(t)  =  r(0)  -  c'Q(t)c;  v(0)=r(0)  (3.1.19) 

r(0)  =  c'Q(0)c  +  +  0^^  =  r(0)  +  (3.1.20) 

Q(t+1)  =  AQ(t)A'  +  v(t)k(t)k'(t);  Q(0)»0  (3.1.21) 

0(0)  =  AQ(0)A'  +  a^j2h(l)h'(l)  (3.1.22) 


Figure  3.1.2  shows  the  Kalman  predictor  with  the  assumed  signal  model, 
the  Innovations  representation  #1,  for  the  noisy  case.  By  inverting 
the  Kalman  predictor  a  new  signal  model  is  obtained.  This  new  model 
is  driven  by  the  innovations  sequence  of  the  Kalman  predictor.  The 
noisy  measurement  process  is  included  in  this  model  since  the  output  is 
z(t)=y(t)+n(t)  not  just  y(t).  Figure  3.1.3  shows  this  new  signal  model 
along  with  the  Kalman  predictor.  Because  of  how  this  model  is  derived 
the  state  vector  of  this  new  model  is  equal  to  the  one  of  the  Kalman 
predictor  is  so  we  can  say  that  the  Kalman  predictor  estimates  the 
state  of  this  new  model  with  no  error  and  reconstructs  the  innovations 
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sequence  perfectly.  But  of  course  the  new  state  vector  is  not  equal 
to  the  state  vector  of  the  old  signal  model. 

3.2  The  Kalman  Filter 

We  said  in  section  3.1  that  (one  step)  prediction  is  the  process 
of  estimating  the  state  vector,  )((t),  of  a  signal  model  given  measure¬ 
ment  data  up  to  time  t-1.  Estimating  the  state  vector  given  data  up 
to  time  t  is  said  to  be  a  filtering  process.  The  Kalman  filter  is  an 
inverse  of  the  Innovations  representation  #2  and  as  in  the  Kalman  pre¬ 
dictor  the  innovations  sequence,  u(t),  plays  the  main  role.  We  con¬ 
sider  two  cases: 

1)  Noise-Free  Measurement  Process:  z(t)  =  y(t) 

2)  Noisy  Measurement  Process:  2(t)  =  y(t)  +  n(t) 

As  before  n(t)  is  the  measurement  noise. 

Noise-free  measurement  processes.  The  signal  model  is  now  the 
Innovations  representation  #2; 

x(t/t)=Ax(t-l/t-l)+fe(t)u(t);  x(-l/-l)=0  (3.2.1) 

y(t-l)=c'x(t-l/t-l)  (3.2.2) 

The  idea  is  to  reconstruct  the  innovations  sequence,  u(t),  from  the 
measurement  data. 

Assume  we  have  a  predicted  output  y(t).  Define  an  innovations 
sequence  as 

u(t)  =  y(t)  -  y(t)  (3.2.3) 

Plug  it  into  equation  (3.2.1): 

x(t/t)  =  Ax(t-l/t-l)  +  fe(t)u(t) 

=  Ax(t-l/t-l)  +  fc(t)[z(t)-y(t)]  (3.2.4) 

Figure  3.2.1  shows  this  process.  We  have  not  said  how  to  get  the 
estimate  y(t).  One  way  is  to  have  the  Kalman  predictor  produce  it  but 
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that  does  not  look  very  attractive  computationally.  First  note  that 
the  state  vector  in  the  Kalman  predictor  can  be  written  as 

2^(t/t-l)  =  Ax^(t-l/t-l)  (3.2.5) 

and  then  the  predicted  output  becomes 

y(t)  =  £‘2^{t/t-l)  =  c*Ax{t-l/t-l)  (3.2.6) 

This  means  we  can  produce  the  estimate  of  the  output  y(t)  using  the 
old  state  >^(t-l/t-l)  in  the  Kalman  filter.  Figure  3.2.2  shows  the 
complete  Kalman  filter.  Like  the  Kalman  predictor  the  Kalman  filter 
reconstructs  the  innovations  sequence,  u(t),  and  estimates  the  state 
of  the  signal  model  (the  Innovations  representation  iZ)  exactly  in  the 
case  of  a  noise-free  measurement  process.  Figure  3.2.3  show  the  Kalman 
filter  along  with  the  assumed  signal  model. 

Note  that  we  get  the  Kalman  predictor  for  free  in  the  Kalman 
filter  since 

>^(t/t-l)  =  A  x(t-l/t-l) 

The  Kalman  filter  is  also  a  whitening  filter  if  we  think  of  u(t)  as  the 


output.  The 

defining  equations  for  the 

Kalman 

filter  in 

the  noise- 

free  case  are 

,  x(t/t) 

=  Ax(t-l/t-l)  +  ^(t)u(t); 

x(-l/- 

1)=0 

(3.2, 

.7) 

u(t)  = 

z(t)  -  y(t) 

(3.2, 

.8) 

y(t)  = 

£'A2<(t-l/t-l ) 

(3.2, 

.9) 

ic(t)v( 

t)  =  [$(0)-A^^(t-l)A']c 

(3.2. 

,10) 

V(t)  = 

r(0)  -  c'A.,Ht-l)A'c;  v(0)  = 

r(0) 

(3.2. 

11) 

Qit)  = 

A4?(t-1)A'  +  v(t)k(t)k'  (t); 

^^(-1) 

=0 

(3.2. 

,12) 

>■(0)  = 

A^(0)A'  +  .’^j2h(0)h' (0) 

(3.2. 

,13) 

r(0)  = 

c'.?(0)c  +  0^^=^ 

(3.2. 

,14) 
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Noisy  measurement  process.  With  additive  noise,  the  input  to 
the  Kalman  filter  is 

z(t)  =  y(t)  +  n(t)  (3,1.15) 

where  n(t)  is  the  measurement  noise.  The  added  noise  only  affects  the 
variance  of  2(t) 

r(k)  =  E[2(t)z' (t+k)]  =  E[(y(t)+n(t))(y(t)+n(t))'] 


=  r(k)  +  o^26(k)  (3.2.16) 

and  for  k=0, 

r(0)  =  r(0)  +  0^2  =  c'5(0)c  +  (3.2.17) 

The  defining  equations  for  the  Kalman  filter  are 

x(t/t)  =  Ax(t-l/t-l)  +£(t)u(t);  x(-l/-l)=0  (3.2.18) 

u(t)  =  z(t)  -  y(t)  (3.2.19) 

y(t)  =  c'Ax(t-l/t-l)  (3.2.20) 

|(t)v(t)  =  [$(0)  -  A$(t-l)A']c  (3.2.21) 

v(t)  =  r(0)  -  c'A$(t-l)A'c;  v(0)*r(0)  (3.2.22) 

$(t)  =  A$(t-1)A'  +  v(t)^(t)k'(t);  $(-l)=0  (3.2.23) 

e(0)  =  A«(0)A'  +  a^j2h(o)h'(o)  (3.2.24) 

?(0)  =  c'e(0)c  +  0^2  (3.2.26) 


Figure,  3.2.4  shows  the  Kalman  filter  with  the  assumed  signal  model  in 
the  noisy  case.  Inverting  the  Kalman  filter  a  new  signal  model  is 
obtained.  This  signal  model  is  driven  by  the  innovations  sequence  of 
the  Kalman  filter  and  includes  the  measurement  noise.  Figure  3.2.5 
shows  this  new  model  along  with  the  Kalman  filter.  The  Kalman  filter 
reconstructs  the  innovations  sequence  of  this  new  model  with  no  error 
and  estimates  the  state  perfectly.  Note  however  that  the  perfectly 
estimated  state  consists  of  imperfect  estimates  of  the  underlying 
sequence  y(t).  Table  3  lists  the  governing  equations  for  the  Kalman 
predictor  and  Kalman  filter. 


-a{n)  .  .  .  .-a(l) 
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3.3  The  Kalman  Gain  Calculations 

The  main  computational  problem  in  implementing  the  Kalman 
predictor  or  filter  is  the  calculation  of  the  Kalman  gain  vectors, 

J<(t)  and  k{t).  The  most  obvious  approach  is  to  use  the  equations 
derived  in  section  3.1  and  3.2.  The  governing  equations  for  the  pre¬ 
dictor  gain,  J<(t)  are 

k(t)v(t)  =  A[0(O)-Q(t)]c'  +  o^2h(l) 
v(t)  =  f(0)  -  c'Q(t)c;  v(0)=r(0) 

Q(t+1)  =  AQ(t)A'  +  v(t)k(t)k'(t);  Q(0)=0 
Q(0)  =  Aq(0)A-  +  a^j2h(l)h'{l) 
r(0)  =  £'Q(0)£  +  ~  >^(0)  +  0^2 

On  the  other  hand  the  governing  equations  for  the  filter  gain,  kit) 
are 

ic(t)v(t)  =  [5(0)-A$(t-l)A']£ 

v(t)  =  r(0)  -  £'A5(t-l)A'£;  v(0)=r(0) 

Qit)  =  A$(t-1)A'  +  v(t)ic(t)fe'(t);  «(-l)=0 
QiO)  =  A«(0)A'  +  o^^2h(0)h’(0) 
r(0)  =  c'QiO)c  +  =  ^(o)  +  0^2 

These  equations  are  of  Ricatti  type  and  are  not  very  attractive  from  a 
computational  point  of  view.  A  simpler  and  faster  algorithm  is  de¬ 
sired. 

Remember  that  the  unit  pulse  response  h^(i)  of  the  two  Inno- 
varions  representations  can  be  written  as 
0  i<0 

h(0)  i=0 

c'A’'^k(t)  i>=l 


h^i)  = 


(3.3.1) 
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or 


>‘(i )  =  I  °  i 

(  c'A^ 


Ht) 


i<0 

i>=0 


(3.3.2) 


Define  the  vector  (i)  as 

h^(i)'  =  [h^(i)  h^(i+l)  .  .  .  h^(i+n-l)]  (3.3.3) 

This  gives  the  following  fundamental  relations: 

k(t)=h^(l)  (3.3.4) 

k{t)  =  h^(0)  (3.3.5) 

All  that  is  necessary  to  get  Jc(t)  and  k{t)  is  to  calculate  the  first 
n+1  of  the  time  varying  impulse  responses  for  each  t.  This  is  shown 
in  Gueguen  and  Scharf  [1]. 

LeRoux  and  Gueguen  [7]  introduced  an  algorithm,  the  so-called 
Impulse  Response  Algorithm,  to  calculate  the  impulse  response  recur¬ 
sively.  By  using  the  fact  that  the  underlying  process  is  an  ARMA(p,q) 
a  substantial  savings  in  calculations  is  possible  as  Dugre'  [8]  points 
out.  This  results  in  the  so-called  Fast  Kalman  Algorithm  of  Morf, 
Sidhu  and  Kailath  [9],  which  they  derived  using  Chandrasekar  type 
equation  for  vector  processes.  The  scalar  version  of  this  fast  algo¬ 
rithm  is  the  same  as  discussed  by  Pearlman  [10]. 

What  follows  is  a  derivation  of  the  Impulse  Response  Algorithm 
for  the  scalar  case.  This  gives  a  fast  way  of  computing  Kalman  gains 
and  leads  to  the  Fast  Kalman  Algorithm.  I  will  follow  the  derivation 
of  Dugre'  [3]  closely  but  extend  it  to  cover  both  k(t)  and  k{t).  He 
derived  the  algorithm  only  for  the  predictor  gain  vector,  j<(t).  Hav¬ 
ing  recursions  for  both  the  predictor  aiid  filter  gain  vectors  bypasses 
the  problem  of  requiring  the  inverse  of  A  to  exist  to  get  the  filter 
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gain  vector  k{t)  as 

k{t)  =  A"^  k(t) 

Impulse  response  algorithm.  The  first  n+1  values  of  the  output 


y(t)  can  be  written  in  matrix  form  as 
Y(t)  =  K'(t)  U(t) 


where 


X'(t)  =  [y(t)  y(t-l) 
U'(t)  =  [u(t)  u(t-l) 


y(i)  y(o)] 

u(l)  u(0)] 


(3.3.6) 


(3.3.7) 

(3.3.8) 


h^-O) 


h^“^(0) 


.  .  .  h'(t-l)  h^(t) 


K‘(t)  = 


(3.3.9) 


h°(l) 

h°(0) 

V.  J 

Then  the  covariance  matrix  R(t)  ®  E[Y(t)Y'(t)]  can  be  written  as  [11] 
R(t)  =  K'(t)  diag[v(t)  .  .  .  v(0)]K(t)  (3.3.10) 


R(t)A(t)  =  K' (t)diag[v(t)  .  .  .  v(0)] 
Here  we  define  A(t)  as  follows: 

’l  0  .  .  . 

a^l)  1 


(3.3.11) 


A(t)=K‘'(t)  = 


(3.3.12) 


[a^(t)  a^’^(t-l)  . 
Now  the  right  hand  term  of  (3.3.11)  is 


a'(l)  1 
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K'  (t)diag[v(t)  .  .  .  v(0]  = 


s‘(0) 


s*(t-l)  s“(t) 


S§(1) 

s“(0) 


(3.3.13) 


where 


s^(i)  =  E[y(t+i)v'(t)]  =  h^(i)v(t) 


(3.3.14) 


The  term  s^(i)  is  a  scaled  version  of  the  impulse  response  h^(i).  It 
is  also  the  cross  correlation  between  the  input  at  t  and  the  output  at 
t+i.  For  i<0  s^(i)  can  be  interpreted  as  an  anti  causal  impulse  re¬ 
sponse.  In  matrix  form  it  satisfies 


r(l+t) 


.  .  r(l) 


r(t+l+t)  . 


r(t+l) 


A(t)  = 


For  example 


s^-t-1) 


s^(-2)  s°(-l) 
s°(-2) 


s^(-k-2)  .  . 


s°(-t-l) 


s*^“^(-k)  =  r(k)+r(k-l)a'^’^  (1)  +  .  .  .  r(l)a'^"^(k-l) 
By  inspecting  equation  (3.3.11)  it  follows  that 

s^(k)  =  r(k)  +  r(k+l)a^(l)  +  ,  .  .  +  r(k+t)a^(t) 


(3.3.15) 


(3.3.16) 


(3.3.17) 
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s^'^^(k)  =  r(k)  +  r(k+l)a^'^^(l)  +  .  .  . 

+  r(kn+l)a^'^^(t+l)  (3.3.18) 

What  we  need  is  a  recursion  for  s^(k). 

Levison-Durbin  algorithm.  The  Levinson-Durbin  algorithm  gives 
a  recursion  for  a^(i)  [12],  [13]: 

a^'^^(i)  =  a^(i)  +  X(t+l)a^  (t+l-i)  i=l,2,  .  .  .,  t  (3.3.19) 

X(t+1)  =  -[r(t+l)+a^(l)r(t)+  .  .  .  +a^(t)r(l)]/v(t)  (3.3.20) 

v(t+l)  =  v(t)[l  -  X2  (t+1)]  (3.3.21) 

a^‘^^(t+l}  =  X(t+1)  and  a^‘^^(0)  =  a^(0)  =  1 
We  can  substitute  for  a^^^(k)  in  equation  (3.3.18)  to  obtain 
s^^^(k)  =  (r(k)+r(k+l)a^(l)+  .  .  .  +r(k+t)a^(t) )  + 

X(t+l)[r(k+l)a^(t)+r(k+2)a^(t-l)+  .  .  .  +r(k+t)a^ 
+r(k+t+q)]  =  s^(k)  +  X(t+l)s^(-t-l-k)  (3.3.22) 

Now  we  need  a  recursion  for  s^(-t-l-k).  It  follows  from  inspection  of 
equation  (3.3.15).  We  write 

s^'^^(-t-l-k)=r(t+l+k)+r(t+k)a^'^^(l)+r(t+k-l)a^'^^(2)  +  .  .  . 
+r(k)a^^^(t+l)  =  [r(t+l+k)+r(t+k)a^(l)+  .  .  . 

+r(k+l)a^(t)]  +  X(t+l)[r(k+t)a^(t)+  .  .  .  +r(k+l)a^(l) 
+r(k)]  =  s^(-t-l-k)  +  X(t+l)s^(k)  (3.3.23) 

Here  we  have  substituted  a^(k)+X(t+l)a^(t+l-k)  for  a^^^(k)  as  before. 
Now  we  can  write  the  recursion  for  s^(k)  as 
s^'^^(k)  =  s^(k)  +  X(t+l)s^(-t-l-k) 
s^'^^(-t-l-k)  =  s^(-t-l-k)  +  X(t+l)s^(k) 

X(t+1)  =  -s^(-t-l)/s^(0) 
with  initial  conditions 

s^(k)  =  r(k)  V  k  and  X(0)  =  -r(l)/r(0) 


(3.3.24) 

(3.3.25) 

(3.3.26) 


This  is  the  so-called  Impulse  Response  Algorithm  [7].  By  substi¬ 
tuting  s^(0)  =  v(t)  and  s^(k)  =  h^(k)v(t),  and  adding  on  the  recursion 
for  a^(k)  from  the  Levinson-Durbin  Algorithm,  we  get  a  simultaneous 
recursion  for  h^(k)  and  a^(k)  as: 

h^'^^(k)  =  [h^(k)  +  X(t+l)h^  (.t-l-k)]/(l  -  X2(t+1))  (3.3.27) 

h^'^^(-t-l-k)=[h^(-t-l-k)+X(t+l)h^(-t-l-k)]/(l-X2(t+l))  (3.3.28) 

X(t+1)  =  -h^(-t-l)v(t)/v(t)  =  -h^(-t-l)  (3.3.29) 

a^'^^(k)  =  a^(k)  +  X(t+l)a^(t+l-k) ,  k=l,2 . t  (3.3.30) 

a^‘^^(t+l)  =  X(t+1)  (3.3.31) 

with  initial  conditions  h^(k)=r(k)/r(0)  and  X(0)=-r(l)/r(0). 

Note  that  we  can  write  s^^^(O)  as 

$^■^^(0)  =  s^(0)(l  -  X2(t+1))  (3.3.32) 

because  of  equation  (3.3.21)  and  s^(0)=v(t). 

An  obvious  drawback  in  using  the  Impulse  Response  Algorithm  is 
that  one  must  know  t  ahead  of  time  to  calculate  s^(-t-l-k)  [8].  This 
can  be  seen  in  the  following  way.  Write  s^(-t-l-k)  as 
s^(-t-l-k)  =  s^‘^(-t-l-k)  +  X(t)s^‘^(k) 

=  s^’^(-t-l-k)  +  X(t-l)s^"^(k)  +  X(t)s^"^(k) 

=s°(-t-l-k)  +  x(t-i)s^'^'''(k) 

i=0 

It  is  clear  that  one  must  start  with  s^(-t-l-k)  to  get  s^(-t-l-k). 

This  depends  on  t.  In  particular  for  k=n  one  must  start  with 
s^(-t-l-n)  to  calculate  s^(-t-l-n)  where  n=max(p,q).  This  drawback 
can  be  corrected  by  noting  the  following: 


[r(n+l+t)  r(n+t)  .  .  .  r(n+l)J  = 


[-a(n) 


-a(2)  -a(l)]  r»-(t+l) 


r(t+n) 


Multiplying  both  sides  by  [1  a^(l)  .  .  .  a^(t)]'  results  in 
s^-t-l-n)  =  -[a(n)  .  .  .  a(l)][s^(-t-l)  .  .  .  s^(-t-n)]'  (3. 

So  there  is  no  need  to  go  all  the  way  down  to  s^(-t-l-n)  to  get 
s^(-t-l-n).  One  can  stop  at  s^(-n-l)  and  proceed  recursively  [8] 
Fast  Kalman  algorithm.  Using  equation  (3.3.24)  we  write 


s‘*'(0) 


s‘(0) 

s‘(l) 


■A(t+1) 


S^(-t-l) 

s^-t-2) 


s^^^(n-l) 


s^(n-l) 


s^(-t-n) 


s^d) 

s^{2) 


-X(t+1) 


s^(-t-2) 

S^(-t-3) 


s^(n) 


s^(-t-n-l) 


In  matrix  form  these  are 

K^^^(O)  =  K^O)  -  (c'LVv(t))L^ 

K^^^(l)  =  K^(l)  -  (c'LVv(t))A 

vhere 

K^O)'  =  [s^(O)s^l)  .  .  .  s^n-1)] 

K^l)'  =  [s^l)s^(2)  .  .  .  s^(n)] 

(lV=  [s^(-t-l)s^(-t-2)  .  .  .  s^(-t-n)] 

(A  L^)'  =  [s^(-t-2)s^(-t-3)  .  .  .  s^(-t-n-l)] 
s^(-t-l)  =  c'L^  and  s^(0)  =  v(t) 


All  we  need  now  is  a  recursion  for  Using  equation  (3.3.25)  we 
write 


•  " 

s^'^^(-t-2) 

s^-t-2) 

s‘(l) 

s^^^(-t-3) 

• 

= 

s^(-t-3) 

-(c'LVv(t)) 

s*(2) 

s^'^^(-t-n-l) 

• 

s^(-t-n-l) 

s‘(n) 

=  AL^  -  (c‘LVv(t))K^(l) 

=  A[L^  -  (c'LVv(t))K^(0)] 

Here  we  have  used  the  fact  that 

K^(l)  =  v(t)k(t)  =  v(t)Afe(t)  =  AK^(O) 

Summary  of  Recursions  for  K^(0)  and  K^(l) 


For  K^(0): 

=  K^(0)  -  (c'LVv(t))L^  (3.3.34) 

=  A  [L^  -  (c'LVv(t))]K^(0)  (3.3.35) 

v(t+l)  =  v(t)[l  -  (c’LVv(t))2]  (3.3.36) 

Initial  conditions; 

K°(0)'  =  [r(0)  r(l)  .  .  .  r(n-l)]  (3.3.37) 

L°(0)'  =  [r(l)  r(2)  .  .  .  r(n)]  (3.3.38) 

v(0)  =  r(0)  (3.3.39) 

For  K^{1): 

jC^+l(l)  =  K^(l)  -  (c'LVv(t))AL^  (3.3.40) 

L^+1  =  AL^  -  (c'L/v{t))K^(l)  (3.3.41) 

v(t+l)  =  v(t)  [1  -  (c'LVv(t))2]  (3.3.42) 

Initial  conditions: 


K°(0)'  =  L°(0)*  =  [r(l)  r(2)  .  .  .  r(n)] 


(3.3.43) 
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The  Fast  Kalman  algorithm  follows  by  substituting  v(t)fe(t)  for 
K^(0),  v(t)k(t)  for  K^(l)  and  v(t)l(t)  for  and  solving  for  k{t) 


and  Jc(t). 

The  Filter  Gain  Vector,  fe(t) 

fe(t+l)  =  Ck(t)  -  (c'i(t))i(t)]/d(t)  (3.3.45) 

l(t+l)  =  A[l(t)  -  (c*l{t))k(t)]/d(t)  (3.3.46) 

d(t)  =  1  -  (c'l(t))2  (3.3.47) 

v(t+l)  =  v(t)d(t)  (3.3.48) 

Initial  conditions  are 

ic'(O)  =  [r(0)  r(l)  .  .  .  r(n-l)]/v(0)  (3.3.49) 

r(0)  =  [r(l)  r(2)  .  .  .  r(n(]/v(0)  (3.3.50) 

v(0)  =  r(0)  (3.3.51) 

The  Predictor  Gain  Vector,  k(t) 

k(t+l)  =  Ck(t)  -  (c*i(t))Al(t)]/d(t)  (3.3.52) 

l(t+l)  =  [Al(t)  -  (c'l(t))k(t)]/d(t)  (3.3.53) 

d(t)  =  1  -  (£'l(t))2  (3.3.54) 

v(t+l)  =  v(t)d(t)  (3.3.55) 

Initial  conditions  are 

k'(0)  =r(0)  =  [r(l)  r(2)  .  .  .  r(n)]/v(0)  (3.3.56) 

v(0)  =  r(0)  (3.3.57) 

In  the  case  of  noisy  measurement  process  then  we  must  use 
v(0)  =  r(0)  =  r(0)  +  (3.3.58) 

instead  of  v(0)=r(0).  Nothing  else  changes. 
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3.4  Properties  of  the  Kalman  Predictor/Filter 

In  section  3.1  and  3.2  we  derived  the  Kalman  predictor  and 
filter  as  inverses  of  the  Innovations  representations  #1  and  #2, 
respectively.  It  should  be  no  surprise  that  the  Kalman  predictor  and 
filter  share  many  of  the  same  properties  as  their  corresponding  signal 
models.  In  this  section  we  will  establish  the  relation  between  the 
state  variance  of  the  Kalman  predictor/filter  and  the  state  variance 
of  the  corresponding  signal  models.  We  also  show  that  the  two  time- 
varying  signal  models  are  asymptotically  time-invariant,  converging  to 
the  Markovian  representations  #1  and  #2,  respectively.  The  Kalman 
predictor  and  filter  converge  to  their  corresponding  time  invariant 
Markovian  predictor  and  filter  structures. 

Consider  first  the  Innovations  representation  #1 
x(t+l/t)  =  Ax(t/t-l)  +  k(t)u(t) 
y(t)  =  c'x(t/t-l)  +  u(t) 
and  its  corresponding  Kalman  predictor 

x(t+l/t)  =  Ax_(t/t-l)  +  £(t)[z{t)-y(t)] 
y(t)  =  c'x(t/t-l) 
z(t)  =  y(t)  +  n(t) 

The  state  variance  of  the  signal  model  Q(t)  satisfies 
q(t+l)  =  A  Q(t)  A'  +  v(t)k(t)k(t);  q(0)=0 
This  can  be  written  as 

q(t)=E[x(t/t-l)x'(t/t-l)]  = 

E[(x(t/t-l)-£(t/t-l)+x(t/t-l))(x(t/t-l)-x(t/t-l)+x(t/t-l))'] 
=E[(x(t/t-l)-x(t/t-l) ) (x(t/t-l )-x(t/t-l) ) • ]+E[x(t/t-l )£• (t/t-1 )] 
=  E(t)  +  q(t)  (3.4.1) 


where  Q{t)  is  the  state  variance  of  the  Kalman  predictor 
Q(t+1)  =  AQ(t)A'  +  v(t)k(t)£'(t);  Q(0)=0 
and  E(t)  is  the  error  variance.  It  follows  that  both  E(t)  and  Q(t) 
are  bounded  by  Q(t)  [2] 

0  <=Q(t)  <=  Q(t)  (3.4.2) 

0  <=E(t)  <=  Q(t)  (3.4.3) 

The  state  variance  of  the  Kalman  predictor  0(t)  is  independent  of  the 
particular  signal  model  used,  depending  only  on  the  covariance  of  the 
measurement  process  z(t).  This  is  so  because  the  Kalman  gain  vector 
J((t)  and  v(t)  depend  only  on  the  covariance  of  the  measurement  process 
r(k)  =  r(k)  +  [2].  On  the  other  hand  the  error  variance  is 

dependent  on  the  particular  signal  model  and  on  the  noise  variance 
of  the  measurement  process.  In  the  case  n(t)=0  the  Kalman  predictor 
estimates  the  state  of  the  Innovations  representation  #1  with  no  error 
and  then 


x(t/t-l)  =  x(t/t-l) 

(3.4.4) 

0 

II 

UJ 

(3.4.5) 

Q(t)  =  Q(t) 

(3.4.6) 

Increasing  the  noise  variance  of  the  measurement  process  results  in  a 
smaller  state  variance  Q(t)  and  larger  error  variance  for  fixed  state 
variance,  Q(t),  of  the  signal  model. 

The  Innovations  representation  #2  is  described  by 
x(t+l/t+l)  =  A2<(t/t)  +  2c(t+l)u(t+l ) 
y(t)  =  c'x(t/t) 

and  its  corresponding  Kalman  filter 

x(t/t)  =  Ax(t-l/t-l)  +  £(t)[z(t)-y(t)] 
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y(t)  =  £'A2((t-l/t-l) 
z(t)  =  y(t)  +  n(t) 

Proceeding  in  a  same  way  as  we  did  before  the  state  variance  of  the 
Innovations  representation  #2  satisfies 

e(t)  =  Ac(t-1)A'  +  v(t)  kit)k'{t)i  C(-1)=0 
This  can  be  written  as 


Q(t)=E[(x(t/t)-£(t/t)+x(t/t))(x(t/t)-x(t/t)+x(t/t))’]  = 

E[ ( X ( t/t ) -X ( t/t ) ) ( X ( t/t )-x ( t/t ) ) ' ]+E[x ( t/t )x ' ( t/t ) ] 

=  £:(t)  +  e(t)  (3.4.7) 

Here  Q(t)  is  the  Kalman  filter  state  variance: 

Q{t)  =  A${t-1)A'  +  v(t)£(t)i'(t);  <?(-l)=0 


and  E{t)  is  the  error  variance.  Both  the  error  and  the  filter  state 
variance  are  bounded  by  the  state  variance  of  the  signal  model  Q(t): 

0  =  E{t)  =  Q{t)  (3.4.8) 

0  =  Q{t)  =  <?(t)  (3.4.9) 

In  case  we  have  n(t)=0  the  Kalman  filter  estimates  the  state  of  the 
Innovations  representation  #2  with  no  error  and  we  have 


x(t/t)  =  x(t/t) 
E{t)  =  0 
Q{t)  =  e(t) 


(3.4.10) 

(3.4.11) 

(3.4.12) 


As  we  mentioned  before  the  two  Innovations  representations  and 
hence  the  Kalman  predictor  and  filter  are  time-varying  but  asymptotic¬ 
ally  time  invariant.  For  the  case  n(t)=0  the  Kalman  predictor  esti¬ 
mates  the  state  of  the  Innovations  representation  #1  perfectly  and  so 
does  the  Kalman  filter  for  the  Innovations  representation  #2.  As  a 
result  we  have 


.  A 


Q(t)  =  Q(t) 

Q{t)  =  e(t) 

Assume  that  the  Kalman  predictor  is  estimating  the  state  of  the 
Markovian  representation  #1  (n(t)=0).  We  have  then 

0(0)  =  >  E(t)  +  Q(t);  E(O)=0(O)  (3.4.13) 

Anderson  and  Moore  [2]  show  that  the  state  variance,  Q(t)  is 
monotone  and  converges  to  a  limit  Q(“)  since  it  is  bounded  by  0(0). 

We  need  to  show  that  this  limit  is  0(0).  If  Q(t)  converges  to  a  limit 
so  does  ^(t)  of  the  Innovations  representation  #2  since  $(t)=AQ(t-l)A' 
Now  the  output  covariance  of  the  Innovations  representation  #2  at 
time  t=<»  is 

r(k)  =  £'A*^$(<»)c 

and  if  it  is  to  match  that  of  the  Markovian  representation  #2  then 
r(k)=£'A*^$(“)c  =  r(k)=£'A*^$(0)£ 
and  it  must  hold  that 

e(«)  =  5(0)  (3.4.14) 

As  a  consequence  since  0(O)=A*5(O)*A'  it  follows  that 

Q(«)  =  0(0)  (3.4.15) 

Also 

lim  v(t)  =  a  2  (3.4.16) 

t^  " 

lim  k(t)  =  h(0)  (3.4.17) 

t-wo 

lim  k{t)  =  ji(l)  (3.4.18) 

t-w 

The  main  point  is  this.  The  Innovations  representations  are  both 
asymptotically  time-invariant  and  converge  to  the  Markovian  represen¬ 
tation  H  and  #2,  respectively.  For  the  case  n(t)=0  the  Kalman 
predictor  and  filter  converge  to  the  Markovian  predictor  and  filter 
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structure.  Figures  3.4.1  and  3.4.2  show  the  Markovian  predictor  and 
filter  structure,  respectively. 

In  the  special  case  the  process  Is  an  AR(p)  instead  of  an 
ARMA(p,q)  then  the  signal  models  converge  in  only  p  steps  (Pearlman 
[10]).  The  same  thing  holds  for  the  Kalman  filter  and  predictor.  The 
reason  for  this  is  that  after  p  steps  the  Impulse  Response  Algorithm 
(see  section  3.3)  has  reached  steady  state: 
hP^l(i)  =  hP(i)  i=0.1,  .  .  ..  p 
v(p+l)  =  v(p) 

since 

A(p+1)  =  -hP(-p-l)  =  aP'^^Cp+l)  =  0 
As  a  consequence 

k(p+l)  =  kip) 
k(p+l)  =  k(p) 
and 


<3(p+l)  =  <3(p) 
Q(p+1)  =  Q(p) 


3.4.2.  The  Markovian  filter. 
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CHAPTER  4 


DIGITAL  IMPLEMENTATION 

In  this  chapter  we  address  some  of  the  problems  encountered  in 
implementing  the  Kalman  predictor  and  Kalman  filter  using  fixed  point 
arithmetic.  The  finite  wordlength  effects  are  (1)  overflow,  (2)  round¬ 
off,  and  (3)  input  and  coefficient  quantization.  A  great  deal  of  work 
has  been  done  to  study  these  effects  and  others  in  time  invariant 
digital  filters.  For  a  review  see  for  example  references  [9]  and 
[10].  Scaling  problems  are  treated  in  [3]  and  [4]  and  scaling  of 
state  space  digital  filters  is  covered  in  [1]  and  [2].  Roundoff  noise 
in  state  space  digital  filters  is  covered  in  references  [1],  [5],  [6], 
and  [7].  Issues  in  digital  implementation  of  control  compensators 
have  been  treated  in  [8]  and  [8].  Many  of  the  results  of  these  studies 
can  be  used  the  Kalman  predictor/filter  problem. 

In  section  4.1  we  look  at  the  scaling  problem.  Section  4.2 
deals  with  the  effects  of  roundoff  on  the  performance  of  the  Kalman 
predictor  and  filter.  Then  in  section  4.3  we  examine  the  computational 
complexity  of  the  predictor/filter  in  terms  number  of  arithmetic 
operations. 

4.1  Scaling  of  Variables 

In  this  section  we  show  how  to  scale  each  element  in  the  state 
vectors  of  the  Kalman  predictor  and  filter  to  limit  the  probability 
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overflow.  The  scaling  rule  used  will  be  the  so-called  1^  scaling  rule 
of  Jackson  [3],  [4]. 

The  set  of  numbers  that  can  be  represented  using  a  fixed  word- 
length  of  m-bits  constitutes  the  dynamic  range  of  variables.  The 

variables  take  on  discrete  values 

ie,  i=0.+/-l.  .  -  (4.1.1) 

where  e  is  the  quantization  step  size.  The  dynamic  range  is  therefore 
bounded  by 

<  ei  <  +e2"'"^  Vi  (4.1.2) 

If  the  bounds  are  forced  to  be  +/-1  then  the  step  size  must  be 

e  =  (4.1.3) 

For  m  =  16  (as  in  a  16-bit  processor)  the  step  size  is 

e  =  2"^^  =  3.05176*10“^  (4.1.4) 

A  filter  is  said  to  be  I2  scaled  if  for  each  variable 

"  Ce2"’‘^]2  (4.1.5) 

Here  f|^(i)  is  the  response  of  the  k^^  variable  to  a  unit  pulse  sequence 
applied  at  the  input,  is  the  variance  of  the  signal  applied  at  the 

input.  0^2  variance  of  the  k^^  filter  state  variable 

when  the  input  is  an  i.i.d.  r.v.  of  variance  The  parameter  6 

determines  the  probability  of  overflow.  Increasing  6  decreases  the 
probability  of  overflow. 

Scaling  of  Variables  in  the  Kalman  Predictor 

The  state  space  equations  for  the  Kalman  predictor  are 
x(t+l/t)  =  Ax(t/t-l)  +  k(t)u(t) 
y(t)  =  c'x(t/t-l) 
u(t)  =  z(t)  -  y(t) 
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i 


Figure  4.1.1  shows  a  block  diagraun  of  the  Kalman  predictor.  To  scale 
the  system  we  introduce  a  diagonal  scaling  matrix  S 


s(l) 


0 


s  = 


(4.1.6) 


and  perform  the  transformation 

x(t/t-l)=S"^  x(t/t-l)  (4.1.7) 

The  scaled  Kalman  predictor  is  then 

x(t+l/t)  =  (S'^  AS)  x(t/t-l)  +  (S“^k(t))u(t)  (4.1.8) 

y(t)  =  (c'S)x(t/t-l)  (4.1.9) 

Figure  4.1.2  shows  the  scaled  Kalman  predictor.  The  Ig  scaling  rule 
for  the  Kalman  predictor  is 

6^  [Q(t)]kk  “  (£2"’"^)'  k»l,2 . n  (4.1.10) 

Here  Q(t)  is  the  state  variance  of  the  Kalman  predictor 

Q(t)  =  v(t)  I  (A’'^k(t-i))(A^'^k(t-i))'  (4.1.11) 

i=l 

These  are  time-varying  scaling  constraints  and  not  very  practical.  A 


better  solution  is  to  find  the  worst  case  scaling  constraints  and  use 
them.  The  state  variance  is  monotonically  increasing  and  converges  in 
the  limit  Q(») 

lim  Q(t)  =  Q(«) 

t-w 

So  the  worst  case  scaling  constraints  are  then 

62  S'^[g(»)]^^  S"''’  =  (e2"’'M2  k=1.2,  .  .  ..  n  (4.1.12) 

which  can  be  written  as 

62  s"^k)q(k,k)  =  (e2'"‘^)2  k=l,2,  . 


.  .  n 


(4.1.13) 


Fig.  4.1.2.  The  scaled  Kalman 


K 


dictor. 
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-2  -2 

where  q(k,k)  =  [Q(<»)]|^,^  and  s  (k)  =  [S  ]|^|^.  The  scale  factors  are 
determined  by 

s(k)  =  k=1.2,  ....n  (4.1.14) 

This  is  the  same  result  obtained  by  Mull  is  and  Roberts  [1].  In  the 
case  (e2'"~^)  =  1  the  scale  factors  are  chosen 

s(k)  =  6/q(k,k)  k=l,2,3,  .  .  .,  n  (4.1.15) 

This  means  the  variance  of  each  scaled  state  element  must  be  less  than 
one  assuming  6  >  1.  A  nice  feature  of  this  scaling  procedure  is  that 
it  allows  us  to  map  the  dynamic  range  of  each  state  element  to  the 
range  of  the  numerical  representation.  The  price  is  an  increased 
number  of  multiplications.  A  simpler  procedure  is  to  choose  all  scale 
factors  equal  to  s 

s  =  max(s(k))  (4.1.16) 

k 

Then  the  scaled  Kalman  predictor  becomes 

x(t+l/t)  =  A)c(t/t-l)  +  (l/s)jc(t)u(t)  (4.1.17) 

y(t)  -  sc'2i(t/t-l)  (4.1.18) 

and  only  (n+1)  multiplications  are  required  for  scaling.  This  is  not 
a  bad  idea  if  all  state  elements  have  variance  of  the  same  order  of 
magnitude.  If  not,  it  is  possible  that  some  state  variables  are  under¬ 
flowing.  By  moving  the  scale  factor  out  in  front  of  the  summation 
node  then  only  one  multiplication  is  needed  for  scaling.  Figure 
4.1.3a  and  4.1.3b  show  the  Kalman  predictor  using  this  simplified 
scaling  procedure. 


Fig.  4.1.3a.  The  scaled  Kalman  predictor.  Simplified  scaling 
procedure. 
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Scaling  of  Variables  in  the  Kalman  Filter 

The  Kalman  filter  is  described  by  the  state  space  equations 
x(t/t)  =  Ax(t-l/t-l)  +  ^(t)u(t) 
y(t)  =  c'Ax(t-l/t-l) 
u(t)  =  2{t)  -  y(t) 

Figure  4.1.4  shows  a  block  diagram  of  the  Kalman  filter.  With  the 
scaling  matrix  S  defined  as  before  the  scaled  Kalman  filter  is 


x(t/t)  =  (S"^AS)x(t-l/t-l)  +  (S‘^fe{t))u(t)  (4.1.19) 

y(t)  =  (£'AS)x(t-l/t-l)  (4.1.20) 

Using  similar  arguments  as  for  the  Kalman  predictor  we  write  down  the 
I2  scaling  constraints  for  the  filter. 

62  s'^(k)q(k.k)  =  (e2"''M2  (4.1.21) 

The  scale  factors  are  then 

s(k)  =  k=l,2.  .  .  .,  n  (4.1.22) 

Here 

q(k.k)  =  (4-1-23) 


where  5(“)  is  the  state  variance  of  the  Kalman  filter  in  steady  state. 

Figure  4.1.5  shows  the  scaled  Kalman  filter.  The  scaling  procedure  is 

simplified  when  all  elements  in  the  matrix  S  are  made  equal  to  s 

s  =  max(s(k))  (4.1.24) 

k 

The  pros  and  cons  of  this  simplification  are  the  same  as  for  the  Kalman 
predictor.  Figure  4.1.6a  shows  the  Kalman  filter  with  this  kind  of 
scaling  and  Figure  4.1.6b  after  having  moved  the  scale  factor  out  in 
front  of  the  summation  node.  Table  4  summarizes  the  basic  scaling 
procedures  derived  in  this  section. 


Fig.  4.1.6b.  The  scaled  Kalman  filter.  Simplified  scaling 
procedure. 
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We  have  seen  how  to  scale  the  state  vectors  in  the  Kalman  predictor 

and  filter.  In  the  simplified  scaling  procedure  we  must  also  make 

* 

sure  that  the  reconstructed  innovations  sequence  u(t)  is  bounded  by 
the  dynamic  range  of  the  numerical  representation.  The  scale  factor 


s  must  be  chosen  larger  or  equal  to 


(4.1.25) 


to  prevent  overflow  in  u(t).  There  is  a  conflict  between  the  scaling 
of  the  state  variables  and  the  scaling  of  the  innovations  sequence 
u(t).  With  increased  measurements  noise  the  state  variance  of  the 
Kalman  predictor/filter  decreases  in  steady  state.  This  means  a 
smaller  scale  factor  s.  However  the  scale  factor  s  must  be  larger  with 
increased  measurement  noise  as  equation  (4.1.25)  shows. 

The  Fast  Kalman  Algorithm  is  very  well  suited  for  fixed  point 
implementation  since  it  is  scaled  by  nature.  This  can  be  seen  in  the 
following  way.  We  know  that  the  gain  vectors  consists  of  the  first 
n+1  values  of  the  time  varying  impulse  response  h^(k).  In  the  deriva¬ 
tion  of  the  Fast  Kalman  Algorithm  we  showed  the  cross  correlation  can 


be  written  as 


s^k)  =  v(t)h^k) 


(4.1.26) 


LeRoux  and  Gueguen  [14]  derived  bounds  on  the  cross  correlations 
ls^(k)l  <=  r(0)  Vk  (4 


Therefore  we  can  write 


(4.1.27) 


(4.1.28 


This  means  all  elements  in  the  Kalman  gain  vector  ^(t)  and  k(t)  are 
bounded  by  +/-1  and  there  is  no  need  to  scale  the  algorithm.  When 
there  is  measurement  noise  the  bound  on  h^(k)  becomes 
lh^k)l  <=  r(0)/(r(0)  +  aj)  <1  H  k 


(4.1.29) 
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4.2  Roundoff  Noise 

The  effects  of  multiplication  roundoff  in  the  Kalman  predictor 
and  Kalman  filter  are  considered  in  this  section.  We  will  be  mainly 
concerned  with  how  the  roundoff  affects  the  state  vectors  since  these 
are  considered  to  be  the  output  of  the  Kalman  predictor  and  filter.  It 
is  also  of  concern  to  know  how  the  roundoff  affects  y(t)  since  this  is 
fed  back  to  the  input  and  any  added  noise  on  y(t)  will  increase  the 
variance  of  u(t).  We  assume  that  the  effect  of  multiplication  roundoff 
can  be  modeled  by  an  additive  white  noise  source  with  zero  mean  and 
variance  z^lVL  and  that  all  such  noise  sources  are  independent  from 
source  to  source.  Here  e  is  the  quantization  step  size  as  before. 

See  for  example  references  [9]  and  [10]. 

Roundoff  Noise  in  the  Kalman  Predictor 

The  state  vector  of  the  Kalman  predictor  is  updated  as 
x(t+l/t)  =  Ax(t/t-l)  +  J<(t)u(t) 

So  for  each  element  in  the  state  vector  x(t+l/t)  there  are  v(i) 
multiplications  required  where 

il  i=l,2,  .  .  . ,  n-1 

(4.2.1) 

n+1  i=n 

This  follows  from  the  structure  of  the  state  matrix  A.  For  every 
multiplication  we  add  one  noise  source  so  the  total  effect  is  an 
additive  white  noise  source  of  zero  mean  and  variance  v(i)e2/i2.  The 
state  update  equation  with  additive  noise  included  is  then 

x(t+l/t)  =  Ax(t/t-l)  +  k(t)u(t)  +  i(t)  (4.2.2) 

where  ^(t)  is  an  uncorrelated  vector  sequence  with  zero  mean  and 
variance  E: 


1 

1 

This  is  the  same  formulation  as  given  by  Mullis  and  Roberts  [1].  The 
state  variance  with  the  effect  of  roundoff  included  is  then 

Q(t)  =  AQ(t-l)A'  +  v(t-l)k(t-l)k'(t-l)  +  E  (4.2.4) 

and  in  steady  state 

Q(<»)  =  AQ(<»)A‘  +  v{<°)k{‘o)k'  (oo)  +  E  (4.2.5) 

The  random  vector  sequence  ^(t)  will  increase  the  state  variance  from 
what  it  was  with  no  roundoff. 

Roundoff  Noise  in  the  Kalman  Filter 

For  the  Kalman  filter  the  state  vector  is  updated  as 
x(t/t)  =  A2<(t-l/t-l)  +  ic(t)u(t) 
which  becomes 

x(t/t)  =  Ax(t-l/t-l)  +  k(t;‘j(t)  +  i(t)  (4.2.6) 

when  the  effect  of  roundoff  is  included.  The  vector  sequence  £(t)  is 
defined  as  before.  The  state  variance  of  the  Kalman  filter  is  then 
Q{t)  =  A5(t-1)A'  +  v(t)k(t)k'(t)  +  E  (4.2.7) 

and  in  steady  state 

$(«)  =  A$(<»)A'  +  v(®)k(®)k' (“)  +  E  (4.2.8) 

The  variance  of  x(t/t-l)  =  Ax(t-l/t-l)  assuming  roundoff  is 
Q(t)  =  Ae(t-1)A' 

=AQ(t-l)A'  +  v(t)k(t)k'(t)  +  AEA'  (4.2.9) 

For  the  Kalman  predictor  assuming  roundoff  we  had 
Q(t)  =  AQ(t-l)A'  +  v(t)k(t)k’(t)  +  E 


This  means  that  using  the  Kalman  filter  as  a  Kalman  predictor  gives 
worse  prediction  results  in  terms  of  state  variance  when  the  effect  of 
roundoff  is  included.  This  increased  state  variance  has  an  effect  on 
the  estimated  output  y(t) 

y(t)  =  £'x(t/t-l)  =  c'[Ax{t-I/t-l)] 
in  terms  of  increased  variance.  In  steady  state  the  variance  of  y(t) 
of  the  Kalman  predictor  is 

E[y(t)y’(t)]  =  c'Q(«>)c 

where  Q(«>)  is  determined  from  equation  (4.2.5).  For  the  Kalman  filter 
the  variance  of  y(t)  in  steady  state  is 
E[y(t)y'(t)]  =  c'qHc 

where  $(«)  is  determined  from  equation  (4.2.9).  This  increased 
variance  of  y(t)  could  affect  the  variance  of  the  innovations  sequence 
u(t)  since  u(t)  =  z(t)-y(t).  To  what  extent  depends  on  the  variance 
of  the  input  sequence  z(t). 

All  derivations  have  been  carried  out  using  the  unsealed  Kalman 
predictor  and  filter.  The  corresponding  results  for  the  scaled  pre¬ 
dictor  and  filter  follow  by  replacing  A  by  (S”^AS),  Jc(t)  by  (S'\(t)), 
fe(t)  by  (S'^?((t))  and  £'  by  (£'S). 

We  have  not  considered  the  effect  of  coefficient  and  input 
quantization.  This  effect  is  basically  deterministic.  We  have  also 
not  determined  how  roundoff  errors  propagate  in  the  calculation  of  the 
Kalman  gain  vectors  using  fixed  point  arithmetic.  It  is  unlikely  that 
the  roundoff  error  can  be  modeled  as  an  additive  white  noise  source 
since  the  gain  values  do  not  change  that  rapidly  between  samples.  In 
Chapter  5  we  give  numerical  results  from  calculating  the  Kalman  gain 
vectors  using  fixed  point  arithmetic. 
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4.3  Computational  Complexity 

This  section  deals  with  the  computational  complexity  of  the 
Kalman  predictor  and  filter  in  terms  of  the  number  of  multiplications, 
divisions,  additions,  and  subtractions.  Memory  requirements  are  also 
considered. 

The  Kalman  predictor  is  described  by  the  following  equations: 
)<(t+l/t)  =  +  k_(t)u(t) 

u(t)  =  z(t)  -  y(t) 
y(t)  =  c'x(t/t-l) 

The  Kalman  filter  is  described  by 

)((t/t)  =  Ax(t-l/t-l)  +  fe(t)u{t) 
u(t)  =  z(t)  -  y(t) 
y(t)  =  c'Ax(t-l/t-l) 

The  computational  burden  can  be  divided  into  the  following: 

1)  Update  the  Kalman  gain  vector  (fe(t)  or  J<(t)) 

2)  Update  the  state  vector  (x(t+l/t)  or  x.(t/t)) 

Even  though  it  is  possible  to  calculate  the  Kalman  gain  vector  ahead 
of  time  we  will  assume  it  is  done  in  real  time  along  with  the  updating 
of  the  state  vector. 

Complexity  of  the  Kalman  Gain  Calculation 
For  the  Kalman  filter  we  have 
k(t+l)  =  [k(t)  -  (c*l(t))l(t)]/d(t) 
i(t+l)  =  A[l(t)  -  (c'l(t))fc(t)]/d(t) 
d(t)  =  1  -  (c'l(t))2 
For  the  Kalman  predictor  we  have 

li(t+l)  =  Ck(t)  -  (c'l(t))Al{t)]/d(t) 


l(t+l)  =  [Al(t)  -  (c'l(t))k(t)]/d(t) 
d(t)  =  1  -  (c'l(t))2 

By  inspection  we  determine  the  number  of  multiplications,  divisions, 
additions,  and  subtractions  required  to  update  the  Kalman  gain  vectors. 
Table  5  gives  the  results.  The  interesting  thing  is  that  it  requires 
a  total  of  (6n+2)  multiplications  and  divisions  to  update  the  predic¬ 
tor  gain  vector  but  only  {5n+2)  to  update  the  filter  gain  vector. 

Using  the  relationship 
k(t)  =  Ak(t) 

to  calculate  the  predictor  gain  vector  requires  the  same  number  of 
arithmetic  operations  as  using  the  defining  equations. 

Complexity  in  Updating  the  State  Vectors 

Again  by  inspection  we  determine  the  number  of  arithmetic 
operations  required  to  update  the  state  vectors.  In  Table  5  we  list 
the  results  along  with  the  total  number  of  arithmetic  operations  needed 
to  go  through  one  cycle  of  the  Kalman  predictor  and  Kalman  filter.  The 
conclusion  is  that  it  takes  a  greater  number  of  arithmetic  operations 
to  implement  the  Kalman  predictor  than  the  Kalman  filter.  The  Kalman 
filter  can  act  both  as  a  filter  and  a  predictor.  So  in  all  respect 
the  Kalman  filter  is  more  efficient  than  the  Kalman  predictor,  except 
that  the  rounding  variance  is  higher  is  the  filter  than  in  the  predic¬ 
tor.  If  the  Kalman  gain  vectot  is  calculated  ahead  of  time  then  the 
Kalman  predictor  and  filter  are  equal  in  terms  of  arithmetic  operations. 
The  number  of  multiplications  and  divisions  required  then  is  on  the 
order  of  (2n). 
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Table  5.  Computational  requirements  for  the  Kalman  predictor  and 
Kalman  filter. 


Updating 


Number  of  Arithmetic  Operations 

#  #  #  # 
Multiplications  Divisions  Additions  Subtractions 


The  Kalman 

Predi ctor 

The  Gain  Vector, 
k(t) 

4n  +  2 

2n 

The  state  Vector, 
x(t+l/t) 

2n 

0 

The  Kalman  Predictor 

6n  +  2 

2n 

The  Kalman  Filter 

The  Gain  Vector, 
kit) 

3n  +  2 

2n 

The  state  Vector, 
x(t/t) 

2n 

0 

The  Kalman  Filter 

5n  +  2 

2n 

Example 


Assume  that  it  takes  ISys  to  multiply  and  divide  and  4us  to  add 
and  subtract.  Then  the  total  cycle  times  for  the  Kalman  predictor  and 
filter  for  the  cases  (1)  n=5  and  (2)  n=10  are  on  the  order  of 

1)  n=5 

Filter:  (7*5+2)*15  +  {l+4)*5  =  580ys 

Predictor:  (8*5+2 )*15  +  (l+4)*5  =  655us 

2)  n=10 

Filter:  (7*10+2)*15  +  {l+4)*5  =  llOSys 

Predictor:  (7*10+2)*15  +  (l+4)*5  =  1255ys 

Memory  Requirements 

For  the  Kalman  predictor  the  storage  requirements  are  on  the 
order  of  (3n+l)  assuming  only  the  state  and  gain  vector  must  be  stored 
along  with  the  denominator  coefficients  a(i)  i=0,n.  When  it  is  re¬ 
quired  to  store  z(t)  and  u(t)  the  total  memory  requirement  is  (3n+3) 
per  time  sample.  The  memory  requirements  for  the  Kalman  filter  are  on 
the  order  of  (4n+3)  assuming  that  the  predictor  state  vector  must  be 
stored  completely. 
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CHAPTER  5 


NUMERICAL  RESULTS 

In  this  chapter  we  give  numerical  results  obtained  by  imple¬ 
menting  the  Kalman  filter  using  floating  point  and  fixed  point 
arithmetic.  The  floating  point  calculations  were  carried  out  on  an 
Intel  Development  System  Model  230  using  Fortran  as  the  program  lan¬ 
guage,  The  fixed  point  implementation  was  carried  out  using  an  Intel 
8086  microprocessor  which  was  a  part  of  the  SDK-86  single  board 
microcomputer.  The  Intel  8086  has  a  wordlength  of  16  bits.  The  actual 
program  was  written  in  assembly  language.  All  the  results  from  the 
floating  point  and  the  fixed  point  implementations  can  be  stored  on  a 
file  for  later  analysis. 

The  input  data  to  the  Kalman  filter  was  generated  by  running 
the  Innovations  representation  #1  with  measurement  noise  added  to  y(t) 
to  obtain  z(t)=y(t)+n(t).  The  data  z(t)  was  generated  on  the  Intel 
Development  System  using  floating  point  arithmetic.  For  information 
on  the  software  and  the  system  configuration  see  Appendices  A  and  B. 

Section  5.1  covers  results  from  the  Kalman  filter  implemented 
using  floating  point  arithmetic.  In  section  5.2  we  look  at  the  cor¬ 
responding  fixed  point  results. 

5. 1  Results  Using  Floating  Point  Arithmetic 

In  this  section  we  present  results  for  a  Kalman  filter  imple¬ 
mented  using  floating  point  arithmetic.  The  example  used  is  an 
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ARMA(3,2)  process: 

_  1  -  1.75  z'^  +  0.8  z~^ _ 

Ti  I?  IT 

1  -  1.5  2  ^  +  1.21  2  ^  -  0.4550  2 

To  start  the  Kalman  filter  we  need  the  denominator  coefficients,  the 
first  4  values  of  the  covariance  sequence  of  the  ARMA  process  and  the 
measurement  noise  variance.  This  is  also  needed  to  generate  the  input 
data  2(t).  The  covariance  sequence  of  the  ARMA  process  is  calculated 
using  the  method  of  Mull  is  and  Roberts  (ref  [5]  in  Chapter  1).  Figure 
5.1.1  shows  a  plot  of  the  first  50  values  of  the  covariance  sequence. 
The  first  4  values  are 

r‘  =  (r(0),r(l),r(2),r(3) 

=  (2.27497,  0.43467,  -1.10293,  -1.14524) 

The  things  we  are  interested  in  are  the  Kalman  filter  gain  vector, 
fe(t),  the  innovations  sequence  u(t)  and  its  calculated  variance  v(t), 
and  the  estimated  output  y(t).  We  also  calculate  the  Kalman  predictor 
gain  vector  £(t)  but  it  is  not  used  for  the  Kalman  filter.  What 
follows  are  results  from  the  Kalman  filter  run  for  three  different 
cases 

1)  0^2  =  0.0,  r(0)/a^2  =  «> 

2)  0^2  =  1.0,  r(0)/o^2  =  2.27497 

3)  0^2  =  2.0,  r(0)/o^2  =  1.13749 

Case  #1 ,  0^2  =  0 

Figure  5.1.2  shows  a  plot  of  the  Kalman  filter  gain  vector 
fe(t)  along  with  the  Kalman  gain  vector  of  the  Innovations  representa¬ 
tion  ff2  ^(t).  The  numbers  from  1  to  3  correspond  to  the  elements  in 
k{X.)  and  the  numbers  4  to  6  correspond  to  the  elements  in  ^(t).  As 
expected  they  are  equal,  i.e.. 
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Fig.  5.1.1.  The  covariance  seq.  r(k)  k  =  0,49  for  an  ARMA(3.2) 
process. 
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fe(t)  =  ic(t) 

In  the  figure  4  overwrites  1,  5  overwrites  2,  .  .  . 

The  initial  value  of  the  Kalman  filter  gain  vector  is 
fe'(O)  =  [r(0)  r(l)  r(2)]/r(0) 

=  (1.00000  0.19107  -0.48481) 

At  time  t=49  the  value  of  the  gain  vector  is 
fe'(49)  =  (1.00000  -0.24994  -0.78495) 

The  final  value  should  be 

k'H  =  (h(0)  h(l)  h(2)) 

=  (1.00000  -0.2500  -0.78500) 

The  h(k)  is  determined  from  equation  (2.2.5). 

Figure  5.1.3  shows  a  plot  of  the  Kalman  predictor  gain  vector 
Jc(t)  along  with  the  Kalman  gain  vector  of  the  Innovations  Representa¬ 
tion  #1  J<(t).  As  can  be  seen  they  are  equal: 

£(t)  =  k(t) 

Again  the  numbers  1,  2  and  3  are  overwritten. 

The  initial  value  of  the  Kalman  predictor  gain  vector  is 
£•(0)  =  [r(l)  r(2)  r(3)]/r(0) 

=  (0.19107  -0.48481  -0.50341) 

At  time  t=49  the  Predictor  gain  vector  is 

£'(49)  =  (-0.24994  -0.78495  -0.42001) 

The  steady  state  value  of  the  Predictor  gain  vector  should  be 
£'(o.)  =  (-0.25000  -0.78500  -0.42000) 

The  anticausal  impulse  response  vector  of  the  signal  model  and 
the  Kalman  filter  is  plotted  in  Figure  5.1.4.  The  numbers  1  to  3 
correspond  to  signal  model  and  the  numbers  4  to  6  to  the  Kalman  filter 
The  two  responses  are  equal.  The  initial  value  is 
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Fig.  5.1.4.  The  anticausal  impulse  response  vectors  Ut)  and  l(t), 


Noise  variance  =  0. 


The  steady  state  value  should  be  zero,  i.e., 
i'(»)  =  (0.0  0.0  0.0) 

At  time  t=49  the  value  is 

r(49)  =  (0.00104  0.00076  -0.00011) 

Figure  5.1.5  shows  a  plot  of  the  output  of  the  signal  model 
y(t)  and  the  estimated  output,  y(t).  The  difference  between  the  two  is 
the  reconstructed  innovations  sequence  since  u(t)=y(t)-y(t) . 

Figure  5.1.6  shows  a  plot  of  the  two  innovations  sequences  of 
the  signal  model  and  the  Kalman  filter,  u(t)  and  u(t).  As  can  be  seen 
the  Kalman  filter  reconstructs  the  innovations  sequence  of  the  signal 
model  exactly.  The  ones  correspond  to  u(t)  and  the  twos  to  u(t).  The 
2s  overwrite  the  Is.  Figure  5.1.7  shows  a  plot  of  the  calculated 
variance  of  u(t)  and  u(t).  The  initial  value  of  v(t)  is 
v(0)  =  r(0)  =  2.27497 
The  final  value  should  be 

v(oc)  =  0^2  =  1.00000 

At  time  t=49  the  calculated  variance  is 
v(49)  =  1.00006 

Case  #2,  =  1.0 

Now  we  have  a  measurement  noise.  Figures  5.1.8  through  5.1.13 
show  the  same  kind  of  information  as  Figures  5.1.2  through  5.1.7.  The 
Kalman  gain  vectors  are  now  different  from  the  noise  free  case.  The 
Kalman  filter  does  not  reconstruct  the  innovations  sequence  of  the 
signal  model  as  before.  The  state  variance  of  the  filter  is  now  less 
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Fig.  5.1.8.  The  Kalman  filter  gain  vectors  k{t)  and  £(t).  Noise 
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Fig.  5.1.13.  The  innovations  variance  v(t)  and  v(t).  Noise 

variance  aj  =  1.0. 
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The  variance  of  the  measurement  noise  is  2.0.  Figures  5.1.14 
through  5.1.19  show  the  same  kind  of  result  as  before.  Table  6  gives 
the  initial  conditions  for  the  variables  we  have  been  looking  at  for 
the  three  different  measurements  noise  level.  The  values  at  time  t=49 
are  also  given. 


5.2  Results  Using  Fixed  Point  Arithmetic 

In  this  section  we  give  results  from  implementing  the  Kalman 
filter  using  fixed  point  arithmetic.  We  compare  the  fixed  point  re¬ 
sults  to  the  corresponding  floating  point  results.  The  same  example  is 
used  as  in  section  5.1: 

H(z)  =  1.  -  1-75  *  0.8  ^ 

1  -  1.5  +  1.21  z"^  -  0.4550  z"'^ 

As  in  section  5.1  we  run  the  fixed  point  version  of  the  Kalman  filter 

for  three  different  cases 


1) 

=  0.0. 

r(0)/aj^2  =  00 

2) 

r(0)/a^2  =  2.27497 

3) 

V  =  2.0. 

r(0)/o^2  =  1.13749 

Kalman 

filter  was 

scaled  using  the  simplified  scaling  procedure 

shown  in  Figure  4.1.6b.  The  purpose  of  the  scaling  is  to  limit  the 
probability  of  overflow  in  the  state  variables.  The  scale  factor  must 
be  chosen  large  enough  to  make  the  scaled  input  z(t)/s  less  than  +/-1. 
The  necessary  scaling  of  the  input  can  be  determined  in  the  same  way  as 
for  the  state  vector.  The  variance  of  z(t)  is  r(0)  +  j^e  scale 

factor  for  the  input  is  chosen  as 


95 


asiMLixi  vmuc 


MIN  VALiiF  >  -.:ase 


MX  MLK  « 


.1  .1 

isnn  ♦ - ♦ — 

1 1 
Z  I 

3  I  3 

4  I  3 

5  I  3 

6  I  3 

7  I  3 
81  3 
»  I  3 
to  13 
11  13 
13  13 
13  13 
M  13 
IS  13 
18  13 

17  13 

18  13 

19  13 

20  13 

21  13 

22  13 

23  !3 

24  13 
3S  13 
28  (3 

27  13 

28  13 

29  13 
31  13 

31  13 

32  13 

33  13 

34  13 

35  13 

38  13 
37  13 
31  13 

39  13 

40  13 

41  13 

42  13 

43  13 

44  13 

45  13 
48  13 

47  13 

48  13 

49  13 
31  13 

INtn  ♦ - ♦ - 

.1  .1 
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In  the  case  =  2.O  the  variance  of  z(t)  becomes  4.27497  and  then 
s  =  2.06766 

Choosing  s  =  4.00  corresponds  to  6  =  1.9  which  should  be  fine. 

The  denominator  coefficients  are  stored  at  their  half  value  since  the 
dynamic  range  is  limited  by  +/*1*0  and  the  largest  denominator  coeffi¬ 
cient  is  -1.5.  The  same  scalefactor  is  used  in  all  three  cases, 
s=4.0.  Also  we  used  truncation  instead  of  rounding. 

Case  #1,  0^2  =  0 

Figures  5.2.1  through  5.2.6  show  plots  of  fe(t),  £(t),  Ut), 
y(t),  u(t),  and  v(t)  calculated  using  floating  and  fixed  point  arith¬ 
metic.  The  lower  numbers  correspond  to  the  floating  point  calculations 
As  can  be  seen  the  fixed  point  results  follow  the  corresponding  float¬ 
ing  point  results  closely.  These  are  very  positive  results. 

Case  #2.  o„2  =  1.0 
n 

Figures  5.2.7  through  5.7.12  show  both  the  floating  and  fixed 
point  results  in  the  case  the  measurement  noise  variance  is  1.0. 

Case  #3,  0^2  =  2.0 

Figures  5.2.13  through  5.2.18  show  the  same  thing  as  Figures 
5.2.7  through  5.2.12.  Now  the  measurement  noise  variance  is  2.0. 

If  anything  the  measurement  noise  seems  to  improve  the  accuracy  of  the 
fixed  point  calculation  compared  with  the  floating  point  ones.  It  is 
also  worth  noting  that  the  rate  of  convergence  increases  with  increased 
measurement  noise.  Compare  for  example  Figures  5.2.1,  5.2.7,  and 
5.2.13.  Table  6  gives  the  initial  values  for  fc(t),  £(t),  X(t),  and 
v(t)  for  both  the  fixed  p  int  and  floating  point  calculations.  The 
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Fig.  5.2.4.  The  estimated  output  y(t)  calculated  using  floating 
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Fig.  5.2.6.  The  innovations  variance  u(t)  calculated  using  floating 
and  fixed  point  arithmetic.  Noise  variance  o^^  =  0. 
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Fig.  5.2.8.  The  Kalman  predictor  gain  vector  £(t)  calculated  using 
floating  and  fixed  point  arithmetic.  Noise  variance  *  l.O. 
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Fig.  5.2.9.  The  anticausal  impulse  response  vector  ^(t)  calculated 
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Fig.  5.2.10.  The  estimated  output  y(t)  calculated  using  floating  and 

fixed  point  arithmetic.  Noise  variance  o  ^  =  l.O. 
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Fig,  5.2.13.  The  Kalman  filter  gain  vector  fc(t)  calculated  using 
floating  and  fixed  point  arithmetic.  Noise  variance  a  ^  =  2.0. 
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Fig.  5.2.14.  The  Kalman  predictor  gain  vector  £(t)  calculated  using 

floating  and  fixed  point  arithmetic.  Noise  variance  o  ^  =  2.0. 
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Fig.  5.2.15.  The  anticausal  impulse  response  vector  l^(t)  calculated 
using  floating  and  fixed  point  arithmetic.  Noise  variance  =  2.0. 
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rig.  5.2.18.  The  innovations  variance  v(t)  calculated  using  floating 

and  fixed  point  arithmetic.  Noise  variance  o  ^  =  2.0. 
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values  at  time  t=49  are  also  given.  One  of  the  major  positive  conclu¬ 
sions  of  this  thesis  is  that  these  fixed  point  results  compare  very 
closely  to  the  floating  point  results.  Measurement  noise  is  desirable 
since  it  enhances  the  performance  of  the  fixed  point  version  of  the 


Kalman  filter. 
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CHAPTER  6 
CONCLUSIONS 

In  this  thesis  we  have  studied  Kalman  filtering  of  ARMA 
processes  both  from  a  theoretical  and  from  a  practical  point  of  view. 
Numerical  results  have  been  presented  for  the  Kalman  filter  imple¬ 
mented  in  floating  point  and  fixed  point  arithmetic. 

The  Kalman  predictor  and  Kalman  filter  are  derived  as  inverses 
to  the  Innovations  representations  #1  and  #2,  respectively.  The 
Innovations  representations  are  linear  time-varying  state  space  models 
which  are  asymptotically  time-invariant.  In  the  limit  the  two  converge 
to  the  so-called  Markovian  representations  #1  and  #2,  respectively. 

As  a  result  the  Kalman  predictor  and  filter  converge  to  the  Markovian 
predictor  and  filter  structure. 

A  very  efficient  recursive  algorithm  called  the  Fast  Kalman 
Algorithm  is  derived  to  calculate  the  Kalman  gain  vector  recursively 
for  both  the  filter  and  the  predictor.  The  number  of  multiplications 
and  divisions  needed  to  update  the  Kalman  filter  gain  vector  is  shown 
to  be  5n+2  and  for  the  predictor  gain  vector  6n+2.  To  update  the  state 
vector  of  the  Kalman  filter  and  predictor  requires  only  2n  multiplica¬ 
tions  so  that  it  is  obvious  that  still  the  Kalman  gain  calculation  is 
the  main  computational  burden  in  implementing  the  Kalman  filter  and 
predictor. 

To  implement  the  Kalman  filter  and  predictor  using  fixed  point 
arithmetic  the  state  variables  must  be  scaled  to  limit  the  probability 
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of  overflow.  A  so-called  I2  scaling  rule  is  used.  For  the  Kalman 
predictor/filter  to  be  properly  scaled  the  variance  of  each  state 
element  must  be  less  than  one  in  magnitude  at  all  time.  Scaling  each 
state  element  properly  requires  a  substantial  number  of  multiplica¬ 
tions.  If  the  variance  of  each  state  element  is  on  the  same  number  of 
magnitude  then  the  scaling  can  be  accomplished  by  only  one  multiplica¬ 
tion.  When  using  this  simplified  scaling  procedure  the  scale  factor  s 
must  be  chosen  large  enough  so  the  scaled  input  z(t)/s  is  less  than 
one  in  magnitude  at  all  time.  This  can  result  in  overscaled  state 
when  there  is  a  large  measurement  noise.  The  Fast  Kalman  Algorithm  is 
scaled  by  nature  and  is  very  well  suited  for  fixed  point  implementa¬ 
tion.  The  numerical  results  confirm  this. 

The  effect  of  multiplication  roundoff  is  also  considered  and 
modeled  by  an  additive  white  noise  vector  in  the  state  update  equation. 

Numerical  results  are  presented  for  implementation  of  the 
Kalman  filter  in  both  floating  point  and  fixed  point  arithmetic.  The 
fixed  point  results  compare  very  well  with  the  floating  point  results 
as  the  figures  in  section  5.2  confirm.  The  fixed  point  results  look 
even  better  considering  that  the  simplified  scaling  procedure  was  used 
and  chopping  instead  of  rounding  was  used.  It  is  also  interesting  to 
note  that  increased  measurement  noise  seems  to  have  a  positive  effect 
on  the  performance  of  the  fixed  point  Kalman  filter. 

The  major  finding  of  this  thesis  is  the  surprisingly  good  per¬ 
formance  of  the  fixed  point  implementation  of  the  Kalman  filter.  The 
Fast  Kalman  Algorithm  is  certainly  efficient  but  is  still  the  major 
computational  burden  in  implementing  the  Kalman  filter. 
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Some  of  the  things  we  have  not  covered  in  this  thesis  but 
should  form  the  basis  of  future  work  are: 

--The  effects  of  roundoff  in  the  Fast  Kalman  Gain  Algorithm. 
--Compute  the  Fast  Kalman  Gain  Algorithm  with  other  Kalman 
gain  algorithms  implemented  in  fixed  point  (Lattice, 
Ricatti ). 

--The  wordlength  issue  in  the  Fast  Kalman  Gain  Algorithm  and 
in  the  Kalman  filter/predictor. 

--Input  and  coefficient  truncation. 

--Limit  cycle  behavior. 
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THE  SOFTWARE 

In  this  appendix  we  describe  the  software  programs  used  in  this 
study  and  describe  how  to  use  them.  There  are  5  main  programs. 

MAINl.  This  program  generates  noisy  data  z(t)=y(t)+n(t).  The 
Innovations  representation  #1  is  used  to  generate  y(t).  The  equations 
used  are 

x(t+l/t)  =  Ax(t/t-l)  +  J((t)u(t) 
x(t+l/t+l)  =  x(t+l/t)  +  fe(t)u(t) 
y(t)  =  c'x(t/t-l) 
k(t)  =  Afe(t) 

The  user  gives  the  names  of  the  files  where  the  following  information 
is  stored: 

i)  the  denominator  coefficients  a{i)  i=0,n 
ii)  the  first  n+1  values  of  the  covariance  sequence  of  the 
ARMA  structure:  r(i)  i=0,n. 

The  user  also  gives  the  length,  nz,  of  the  data  sequence  to  be  gen¬ 
erated  and  the  noise  variance,  Varn,  of  the  measurement  noise  n(t). 

At  the  end  there  is  an  option  of  saving  all  the  data  on  a 
single  file  or  all  or  part  of  the  data  on  separate  files.  The  data  to 
be  saved  includes  z(t),  y(t),  u(t),  v(t),  ^(t),  k(t),  l(t),  x(t/t-l), 
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MAIN2.  This  program  performs  floating  point  Kalman  filtering 
of  the  data  sequence  z(t).  The  equations  used  are 
x(t/t)  =  x(t/t-l)  +  fe(t)(z(t)-y(t)) 
x(t/t-l)  =  Ax(t-l/t-l) 
y(t)  =  c'x(t/t-l) 
k(t)  =  ^k{t) 

The  user  gives  the  names  of  the  files  where  the  following  information 
is  stored: 

i)  the  denominator  coefficients  a(i)  i=0,n  , 

ii)  the  first  n+1  values  of  the  covariance  sequence  of  the  ARMA 
structure:  r(i)  i=0,n  , 

iii)  the  noisy  data  sequence  z(t). 

The  user  also  specifies  how  many  values  of  the  noisy  sequence  are 
desired  and  gives  the  assumed  noise  variance,  Varn,  of  the  measurement 
noise  n(t). 

At  the  end  there  is  an  option  of  saving  all  the  data  on  a  single 
file  or  all  or  part  of  the  data  on  separate  files.  The  data  saved  can 
include  z(t),  y(t),  u(t),  v(t),  fc(t),  R(t),  l(t),  x(t/t-l),  and 
x(t/t). 

MAIN3.  This  program  communicates  with  the  MAIN86  which  imple¬ 
ments  the  Kalman  filter  using  fixed  point  arithmetic  on  Intel's  SDK86 
single  board  computer.  Before  running  the  MAIN3  program  on  the  Devel¬ 
opment  system  the  hex  version  of  MAIN86,  MAIN86.HEX  should  be  down¬ 
loaded  first  and  started. 

The  user  gives  the  names  of  the  files  where  the  following 
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i)  the  denominator  coefficients  a(i)  i=0,n  , 

ii)  the  first  n+1  values  of  the  covariance  sequence  of  the 
ARMA  structure:  r(i)  i=0,n  , 

iii)  the  noisy  data  z(t). 

The  user  also  specifies  how  many  values  of  the  2(t}  sequence  are 
desired  and  gives  the  assumed  noise  variance,  Varn,  of  the  measurement 
noise  n(t).  The  scale  factors  Scale_A  and  Scale_Z  for  the  denominator 
coefficients  and  the  noisy  data  sequence  must  also  be  given.  The 
scale  factor  Scale_A  should  be  chosen  such  that  all  the  scaled  denom¬ 
inator  coefficients  will  be  less  than  one  in  magnitude.  The  scale 
factor  Scale_Z  should  be  chosen  such  that  the  scaled  data  sequence 
(Scale_Z)z(t)  and  all  state  elements  are  less  than  one  in  magnitude  at 
all  times. 

At  the  end  the  user  has  the  option  of  saving  all  the  data  on  a 
single  file  or  all  or  part  of  the  data  on  separate  files.  The  data 
saved  can  include  z(t),  y(t),  u(t),  v(t),  C(t),  £(t),  l{t),  x(t/t-l), 
and  x(t/t). 

MAIN4.  This  program  allows  the  user  to  create  a  data  file  and 
generate  the  covariance  sequence  r{i)  i=0,nr  corresponding  to  an  ARMA 
process  when  given  the  file  names  where  the  denominator  and  nominator 
coefficients,  a(i)  i=0,n  and  b(i)  i=0,n,  are  stored.  The  covariance 
is  calculated  using  a  method  of  Mullis  and  Roberts. 

MAIN86.  This  is  the  fixed  point  implementation  of  the  Kalman 
filter.  The  program  is  written  as  assembly  language  and  runs  on  the 
Intel  single  board  computer  S0K86.  See  also  the  description  of  the 
MAIN3  program  which  communicates  with  the  MAIN86.  The  hex  version  of 
MAIN86  must  be  downloaded  to  the  SDK86  from  the  Development  system  and 


executed  before  running  the  fiAIN3  program.  All  the  data  generated  by 
this  MAIN86  program  is  sent  to  the  development  system  so  that  it  can 
be  saved  on  a  file(s)  for  later  analysis.  To  save  on  the  number  of 
arithmetic  operations  we  use  truncation  instead  of  rounding. 
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SYSTEM  CONFIGURATION 

The  system  configuration  under  which  the  software  programs 
described  in  Appendix  A  run  is  the  following: 

1.  Intellect  Series  II  Model  230  Microcomputer  Development 
System  (MDS-230)  which  includes: 

i )  64  k  bytes  RAM 

ii)  Integral  CRT  and  detachable  keyboard 

iii)  Two  flexible  disk  drives  (1  million  bytes) 

2.  SDK-86  single  board  microcomputer  system  with  necessary 
software  and  hardware  to  interface  to  the  development 
system  (SDK-C86). 

3.  A  line  printer. 

Figure  1  shows  the  basic  configuration.  For  more  information  see  for 
example  the  following  reference  manuals: 

A  Guide  to  Microcomputer  Development  Systems. 

Microcomputer  Development  Packages. 

The  8086  Family  User's  Manual  (Appendix  B,  f'evelopment  Tools). 
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